{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Logistic Regression with CmdStanPy\n",
    "\n",
    "TODO: Work in progress \n",
    "\n",
    "Authors: Jonah Gabry, Ben Goodrich, Aki Vehtari, Tuomas Sivula\n",
    "\n",
    "The introduction to Bayesian logistic regression is from a [CRAN vignette](https://cran.r-project.org/web/packages/rstanarm/vignettes/binomial.html) by Jonah Gabry and Ben Goodrich. CRAN vignette was modified to a [R notebook](https://github.com/avehtari/BDA_R_demos/blob/master/demos_rstan/diabetes.Rmd) by Aki Vehtari.  Instead of wells data in CRAN vignette, Pima Indians data is used. The end of the notebook differs significantly from the CRAN vignette.  The R notebook was ported to this Python notebook by Aki Vehtari and Tuomas Sivula."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "This vignette explains how to estimate generalized linear models (GLMs) for binary (Bernoulli) response variables using CmdStanPy.\n",
    "\n",
    "The four steps of a Bayesian analysis are\n",
    "\n",
    "1. Specify a joint distribution for the outcome(s) and all the unknowns, which typically takes the form of a marginal prior distribution for the unknowns multiplied by a likelihood for the outcome(s) conditional on the unknowns. This joint distribution is proportional to a posterior distribution of the unknowns conditional on the observed data\n",
    "2. Draw from posterior distribution using Markov Chain Monte Carlo (MCMC).\n",
    "3. Evaluate how well the model fits the data and possibly revise the model.\n",
    "4. Draw from the posterior predictive distribution of the outcome(s) given interesting values of the predictors in order to visualize how a manipulation of a predictor affects (a function of) the outcome(s).\n",
    "This notebook demonstrates Steps 1-3 when the likelihood is the product of conditionally independent binomial distributions (possibly with only one trial per observation).\n",
    "\n",
    "### Likelihood\n",
    "\n",
    "For a binomial GLM the likelihood for one observation $y$ can be written as a conditionally binomial PMF $$\\binom{n}{y} \\pi^{y} (1 - \\pi)^{n - y},$$ where $n$ is the known number of trials, $\\pi = g^{-1}(\\eta)$ is the probability of success and $\\eta = \\alpha + \\mathbf{x}^\\top \\boldsymbol{\\beta}$ is a linear predictor. For a sample of size $N$, the likelihood of the entire sample is the product of $N$ individual likelihood contributions.\n",
    "\n",
    "Because $\\pi$ is a probability, for a binomial model the link function $g$ maps between the unit interval (the support of $\\pi$) and the set of all real numbers $\\mathbb{R}$. When applied to a linear predictor $\\eta$ with values in $\\mathbb{R}$, the inverse link function $g^{-1}(\\eta)$ therefore returns a valid probability between 0 and 1.\n",
    "\n",
    "The two most common link functions used for binomial GLMs are the logit and probit functions. With the logit (or log-odds) link function $g(x) = \\ln{\\left(\\frac{x}{1-x}\\right)}$, the likelihood for a single observation becomes\n",
    "\n",
    "$$\\binom{n}{y}\\left(\\text{logit}^{-1}(\\eta)\\right)^y \\left(1 - \\text{logit}^{-1}(\\eta)\\right)^{n-y} = \\binom{n}{y} \\left(\\frac{e^{\\eta}}{1 + e^{\\eta}}\\right)^{y} \\left(\\frac{1}{1 + e^{\\eta}}\\right)^{n - y}$$\n",
    "\n",
    "and the probit link function $g(x) = \\Phi^{-1}(x)$ yields the likelihood\n",
    "\n",
    "$$\\binom{n}{y} \\left(\\Phi(\\eta)\\right)^{y} \\left(1 - \\Phi(\\eta)\\right)^{n - y},$$\n",
    "\n",
    "where $\\Phi$ is the CDF of the standard normal distribution. The differences between the logit and probit functions are minor and -- if, as rstanarm does by default, the probit is scaled so its slope at the origin matches the logit's -- the two link functions should yield similar results. Unless the user has a specific reason to prefer the probit link, we recommend the logit simply because it will be slightly faster and more numerically stable.\n",
    "\n",
    "In theory, there are infinitely many possible link functions, although in practice only a few are typically used. \n",
    "\n",
    "\n",
    "### Priors\n",
    "\n",
    "A full Bayesian analysis requires specifying prior distributions $f(\\alpha)$ and $f(\\boldsymbol{\\beta})$ for the intercept and vector of regression coefficients. \n",
    "\n",
    "As an example, suppose we have $K$ predictors and believe --- prior to seeing the data --- that $\\alpha, \\beta_1, \\dots, \\beta_K$ are as likely to be positive as they are to be negative, but are highly unlikely to be far from zero. These beliefs can be represented by normal distributions with mean zero and a small scale (standard deviation).\n",
    "\n",
    "If, on the other hand, we have less a priori confidence that the parameters will be close to zero then we could use a larger scale for the normal distribution and/or a distribution with heavier tails than the normal like the Student's $t$ distribution.\n",
    "\n",
    "### Posterior\n",
    "\n",
    "With independent prior distributions, the joint posterior distribution for $\\alpha$ and $\\boldsymbol{\\beta}$ is proportional to the product of the priors and the $N$ likelihood contributions:\n",
    "\n",
    "$$f\\left(\\alpha,\\boldsymbol{\\beta} | \\mathbf{y},\\mathbf{X}\\right) \\propto f\\left(\\alpha\\right) \\times \\prod_{k=1}^K f\\left(\\beta_k\\right) \\times \\prod_{i=1}^N { g^{-1}\\left(\\eta_i\\right)^{y_i} \\left(1 - g^{-1}\\left(\\eta_i\\right)\\right)^{n_i-y_i}}.$$\n",
    "\n",
    "This is posterior distribution that CmdStanPy will draw from when using MCMC."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logistic Regression Example\n",
    "\n",
    "When the logit link function is used the model is often referred to as a logistic regression model (the inverse logit function is the CDF of the standard logistic distribution). As an example, here we will show how to carry out a analysis for Pima Indians data set similar to analysis from Chapter 5.4 of Gelman and Hill (2007) using CmdStanPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os, sys\n",
    "\n",
    "import arviz as az\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# import stan interface\n",
    "from cmdstanpy import CmdStanModel\n",
    "\n",
    "# add utilities directory to path\n",
    "# so we can import check_utility\n",
    "util_path = os.path.abspath(os.path.join(os.path.pardir, 'utilities_and_data'))\n",
    "if util_path not in sys.path and os.path.exists(util_path):\n",
    "    sys.path.insert(0, util_path)\n",
    "from check_utility import check_div, check_energy, check_treedepth\n",
    "\n",
    "az.style.use('arviz-white')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data\n",
    "\n",
    "First we load and pre-process data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 768 entries, 0 to 767\n",
      "Data columns (total 9 columns):\n",
      " #   Column                    Non-Null Count  Dtype  \n",
      "---  ------                    --------------  -----  \n",
      " 0   Pregnancies               768 non-null    int64  \n",
      " 1   Glucose                   768 non-null    int64  \n",
      " 2   BloodPressure             768 non-null    int64  \n",
      " 3   SkinThickness             768 non-null    int64  \n",
      " 4   Insulin                   768 non-null    int64  \n",
      " 5   BMI                       768 non-null    float64\n",
      " 6   DiabetesPedigreeFunction  768 non-null    float64\n",
      " 7   Age                       768 non-null    int64  \n",
      " 8   Outcome                   768 non-null    int64  \n",
      "dtypes: float64(2), int64(7)\n",
      "memory usage: 54.1 KB\n"
     ]
    }
   ],
   "source": [
    "# load data\n",
    "data_path = os.path.abspath(\n",
    "    os.path.join(\n",
    "        os.path.pardir,\n",
    "        'utilities_and_data',\n",
    "        'diabetes.csv'\n",
    "    )\n",
    ")\n",
    "data = pd.read_csv(data_path)\n",
    "# print some basic info()\n",
    "data.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Pregnancies</th>\n",
       "      <th>Glucose</th>\n",
       "      <th>BloodPressure</th>\n",
       "      <th>SkinThickness</th>\n",
       "      <th>Insulin</th>\n",
       "      <th>BMI</th>\n",
       "      <th>DiabetesPedigreeFunction</th>\n",
       "      <th>Age</th>\n",
       "      <th>Outcome</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>6</td>\n",
       "      <td>148</td>\n",
       "      <td>72</td>\n",
       "      <td>35</td>\n",
       "      <td>0</td>\n",
       "      <td>33.6</td>\n",
       "      <td>0.627</td>\n",
       "      <td>50</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>85</td>\n",
       "      <td>66</td>\n",
       "      <td>29</td>\n",
       "      <td>0</td>\n",
       "      <td>26.6</td>\n",
       "      <td>0.351</td>\n",
       "      <td>31</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>8</td>\n",
       "      <td>183</td>\n",
       "      <td>64</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>23.3</td>\n",
       "      <td>0.672</td>\n",
       "      <td>32</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>89</td>\n",
       "      <td>66</td>\n",
       "      <td>23</td>\n",
       "      <td>94</td>\n",
       "      <td>28.1</td>\n",
       "      <td>0.167</td>\n",
       "      <td>21</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>137</td>\n",
       "      <td>40</td>\n",
       "      <td>35</td>\n",
       "      <td>168</td>\n",
       "      <td>43.1</td>\n",
       "      <td>2.288</td>\n",
       "      <td>33</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Pregnancies  Glucose  BloodPressure  SkinThickness  Insulin   BMI  \\\n",
       "0            6      148             72             35        0  33.6   \n",
       "1            1       85             66             29        0  26.6   \n",
       "2            8      183             64              0        0  23.3   \n",
       "3            1       89             66             23       94  28.1   \n",
       "4            0      137             40             35      168  43.1   \n",
       "\n",
       "   DiabetesPedigreeFunction  Age  Outcome  \n",
       "0                     0.627   50        1  \n",
       "1                     0.351   31        0  \n",
       "2                     0.672   32        1  \n",
       "3                     0.167   21        0  \n",
       "4                     2.288   33        1  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# preview some first rows\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Pregnancies</th>\n",
       "      <th>Glucose</th>\n",
       "      <th>BloodPressure</th>\n",
       "      <th>SkinThickness</th>\n",
       "      <th>Insulin</th>\n",
       "      <th>BMI</th>\n",
       "      <th>DiabetesPedigreeFunction</th>\n",
       "      <th>Age</th>\n",
       "      <th>Outcome</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>768.000000</td>\n",
       "      <td>768.000000</td>\n",
       "      <td>768.000000</td>\n",
       "      <td>768.000000</td>\n",
       "      <td>768.000000</td>\n",
       "      <td>768.000000</td>\n",
       "      <td>768.000000</td>\n",
       "      <td>768.000000</td>\n",
       "      <td>768.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>3.845052</td>\n",
       "      <td>120.894531</td>\n",
       "      <td>69.105469</td>\n",
       "      <td>20.536458</td>\n",
       "      <td>79.799479</td>\n",
       "      <td>31.992578</td>\n",
       "      <td>0.471876</td>\n",
       "      <td>33.240885</td>\n",
       "      <td>0.348958</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>3.369578</td>\n",
       "      <td>31.972618</td>\n",
       "      <td>19.355807</td>\n",
       "      <td>15.952218</td>\n",
       "      <td>115.244002</td>\n",
       "      <td>7.884160</td>\n",
       "      <td>0.331329</td>\n",
       "      <td>11.760232</td>\n",
       "      <td>0.476951</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.078000</td>\n",
       "      <td>21.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>99.000000</td>\n",
       "      <td>62.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>27.300000</td>\n",
       "      <td>0.243750</td>\n",
       "      <td>24.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>3.000000</td>\n",
       "      <td>117.000000</td>\n",
       "      <td>72.000000</td>\n",
       "      <td>23.000000</td>\n",
       "      <td>30.500000</td>\n",
       "      <td>32.000000</td>\n",
       "      <td>0.372500</td>\n",
       "      <td>29.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>6.000000</td>\n",
       "      <td>140.250000</td>\n",
       "      <td>80.000000</td>\n",
       "      <td>32.000000</td>\n",
       "      <td>127.250000</td>\n",
       "      <td>36.600000</td>\n",
       "      <td>0.626250</td>\n",
       "      <td>41.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>17.000000</td>\n",
       "      <td>199.000000</td>\n",
       "      <td>122.000000</td>\n",
       "      <td>99.000000</td>\n",
       "      <td>846.000000</td>\n",
       "      <td>67.100000</td>\n",
       "      <td>2.420000</td>\n",
       "      <td>81.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Pregnancies     Glucose  BloodPressure  SkinThickness     Insulin  \\\n",
       "count   768.000000  768.000000     768.000000     768.000000  768.000000   \n",
       "mean      3.845052  120.894531      69.105469      20.536458   79.799479   \n",
       "std       3.369578   31.972618      19.355807      15.952218  115.244002   \n",
       "min       0.000000    0.000000       0.000000       0.000000    0.000000   \n",
       "25%       1.000000   99.000000      62.000000       0.000000    0.000000   \n",
       "50%       3.000000  117.000000      72.000000      23.000000   30.500000   \n",
       "75%       6.000000  140.250000      80.000000      32.000000  127.250000   \n",
       "max      17.000000  199.000000     122.000000      99.000000  846.000000   \n",
       "\n",
       "              BMI  DiabetesPedigreeFunction         Age     Outcome  \n",
       "count  768.000000                768.000000  768.000000  768.000000  \n",
       "mean    31.992578                  0.471876   33.240885    0.348958  \n",
       "std      7.884160                  0.331329   11.760232    0.476951  \n",
       "min      0.000000                  0.078000   21.000000    0.000000  \n",
       "25%     27.300000                  0.243750   24.000000    0.000000  \n",
       "50%     32.000000                  0.372500   29.000000    0.000000  \n",
       "75%     36.600000                  0.626250   41.000000    1.000000  \n",
       "max     67.100000                  2.420000   81.000000    1.000000  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# some summary\n",
    "data.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preprocess data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# modify the data column names slightly for easier typing\n",
    "# rename DiabetesPedigreeFunction to dpf\n",
    "data.rename(columns={'DiabetesPedigreeFunction': 'dpf'}, inplace=True)\n",
    "# make lower\n",
    "data.rename(columns=lambda old_name: old_name.lower(), inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# removing those observation rows with 0 in selected variables\n",
    "normed_predictors = [\n",
    "    'glucose',\n",
    "    'bloodpressure',\n",
    "    'skinthickness',\n",
    "    'insulin',\n",
    "    'bmi'\n",
    "]\n",
    "data = data[(data[normed_predictors] != 0).all(axis=1)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_15399/3539175312.py:3: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '3     -2.30102\n",
      "4     -3.30102\n",
      "6     -0.30102\n",
      "8     -1.30102\n",
      "13    -2.30102\n",
      "        ...   \n",
      "753   -3.30102\n",
      "755   -2.30102\n",
      "760   -1.30102\n",
      "763    6.69898\n",
      "765    1.69898\n",
      "Name: pregnancies, Length: 392, dtype: float64' has dtype incompatible with int64, please explicitly cast to a compatible dtype first.\n",
      "  data.iloc[:,:-1] -= data.iloc[:,:-1].mean()\n",
      "/tmp/ipykernel_15399/3539175312.py:3: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '3     -33.627551\n",
      "4      14.372449\n",
      "6     -44.627551\n",
      "8      74.372449\n",
      "13     66.372449\n",
      "         ...    \n",
      "753    58.372449\n",
      "755     5.372449\n",
      "760   -34.627551\n",
      "763   -21.627551\n",
      "765    -1.627551\n",
      "Name: glucose, Length: 392, dtype: float64' has dtype incompatible with int64, please explicitly cast to a compatible dtype first.\n",
      "  data.iloc[:,:-1] -= data.iloc[:,:-1].mean()\n",
      "/tmp/ipykernel_15399/3539175312.py:3: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '3      -4.663265\n",
      "4     -30.663265\n",
      "6     -20.663265\n",
      "8      -0.663265\n",
      "13    -10.663265\n",
      "         ...    \n",
      "753    17.336735\n",
      "755    17.336735\n",
      "760   -12.663265\n",
      "763     5.336735\n",
      "765     1.336735\n",
      "Name: bloodpressure, Length: 392, dtype: float64' has dtype incompatible with int64, please explicitly cast to a compatible dtype first.\n",
      "  data.iloc[:,:-1] -= data.iloc[:,:-1].mean()\n",
      "/tmp/ipykernel_15399/3539175312.py:3: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '3      -6.145408\n",
      "4       5.854592\n",
      "6       2.854592\n",
      "8      15.854592\n",
      "13     -6.145408\n",
      "         ...    \n",
      "753    14.854592\n",
      "755     9.854592\n",
      "760    -3.145408\n",
      "763    18.854592\n",
      "765    -6.145408\n",
      "Name: skinthickness, Length: 392, dtype: float64' has dtype incompatible with int64, please explicitly cast to a compatible dtype first.\n",
      "  data.iloc[:,:-1] -= data.iloc[:,:-1].mean()\n",
      "/tmp/ipykernel_15399/3539175312.py:3: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '3      -62.056122\n",
      "4       11.943878\n",
      "6      -68.056122\n",
      "8      386.943878\n",
      "13     689.943878\n",
      "          ...    \n",
      "753    353.943878\n",
      "755    -46.056122\n",
      "760   -140.056122\n",
      "763     23.943878\n",
      "765    -44.056122\n",
      "Name: insulin, Length: 392, dtype: float64' has dtype incompatible with int64, please explicitly cast to a compatible dtype first.\n",
      "  data.iloc[:,:-1] -= data.iloc[:,:-1].mean()\n",
      "/tmp/ipykernel_15399/3539175312.py:3: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '3      -9.864796\n",
      "4       2.135204\n",
      "6      -4.864796\n",
      "8      22.135204\n",
      "13     28.135204\n",
      "         ...    \n",
      "753    -4.864796\n",
      "755     6.135204\n",
      "760    -8.864796\n",
      "763    32.135204\n",
      "765    -0.864796\n",
      "Name: age, Length: 392, dtype: float64' has dtype incompatible with int64, please explicitly cast to a compatible dtype first.\n",
      "  data.iloc[:,:-1] -= data.iloc[:,:-1].mean()\n"
     ]
    }
   ],
   "source": [
    "# scale the covariates for easier comparison of coefficient posteriors\n",
    "# N.B. int columns turn into floats\n",
    "data.iloc[:,:-1] -= data.iloc[:,:-1].mean()\n",
    "data.iloc[:,:-1] /= data.iloc[:,:-1].std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pregnancies</th>\n",
       "      <th>glucose</th>\n",
       "      <th>bloodpressure</th>\n",
       "      <th>skinthickness</th>\n",
       "      <th>insulin</th>\n",
       "      <th>bmi</th>\n",
       "      <th>dpf</th>\n",
       "      <th>age</th>\n",
       "      <th>outcome</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.716511</td>\n",
       "      <td>-1.089653</td>\n",
       "      <td>-0.373178</td>\n",
       "      <td>-0.584363</td>\n",
       "      <td>-0.522175</td>\n",
       "      <td>-0.709514</td>\n",
       "      <td>-1.030559</td>\n",
       "      <td>-0.967063</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-1.027899</td>\n",
       "      <td>0.465719</td>\n",
       "      <td>-2.453828</td>\n",
       "      <td>0.556709</td>\n",
       "      <td>0.100502</td>\n",
       "      <td>1.424909</td>\n",
       "      <td>5.108582</td>\n",
       "      <td>0.209318</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>-0.093734</td>\n",
       "      <td>-1.446093</td>\n",
       "      <td>-1.653578</td>\n",
       "      <td>0.271441</td>\n",
       "      <td>-0.572662</td>\n",
       "      <td>-0.296859</td>\n",
       "      <td>-0.796108</td>\n",
       "      <td>-0.476904</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>-0.405123</td>\n",
       "      <td>2.409934</td>\n",
       "      <td>-0.053078</td>\n",
       "      <td>1.507603</td>\n",
       "      <td>3.255961</td>\n",
       "      <td>-0.368007</td>\n",
       "      <td>-1.056609</td>\n",
       "      <td>2.169953</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>-0.716511</td>\n",
       "      <td>2.150705</td>\n",
       "      <td>-0.853328</td>\n",
       "      <td>-0.584363</td>\n",
       "      <td>5.805571</td>\n",
       "      <td>-0.424924</td>\n",
       "      <td>-0.361940</td>\n",
       "      <td>2.758143</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    pregnancies   glucose  bloodpressure  skinthickness   insulin       bmi  \\\n",
       "3     -0.716511 -1.089653      -0.373178      -0.584363 -0.522175 -0.709514   \n",
       "4     -1.027899  0.465719      -2.453828       0.556709  0.100502  1.424909   \n",
       "6     -0.093734 -1.446093      -1.653578       0.271441 -0.572662 -0.296859   \n",
       "8     -0.405123  2.409934      -0.053078       1.507603  3.255961 -0.368007   \n",
       "13    -0.716511  2.150705      -0.853328      -0.584363  5.805571 -0.424924   \n",
       "\n",
       "         dpf       age  outcome  \n",
       "3  -1.030559 -0.967063        0  \n",
       "4   5.108582  0.209318        1  \n",
       "6  -0.796108 -0.476904        1  \n",
       "8  -1.056609  2.169953        1  \n",
       "13 -0.361940  2.758143        1  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# preview some first rows againg\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# preparing the inputs\n",
    "X = data.iloc[:,:-1]\n",
    "y = data.iloc[:,-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of observations = 392\n",
      "number of predictors = 8\n"
     ]
    }
   ],
   "source": [
    "# get shape into variables\n",
    "n, p = X.shape\n",
    "print('number of observations = {}'.format(n))\n",
    "print('number of predictors = {}'.format(p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stan model code for logistic regression\n",
    "\n",
    "Logistic regression with Student's $t$ prior as discussed above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/**\n",
      " * Logistic regression t-prior\n",
      " *\n",
      " * Priors:\n",
      " *     weights - student t\n",
      " *     intercept - student t\n",
      " */\n",
      "data {\n",
      "    int<lower=0> n;                      // number of data points\n",
      "    int<lower=1> d;                      // explanatory variable dimension\n",
      "    matrix[n, d] X;                      // explanatory variable\n",
      "    array[n] int<lower=0, upper=1> y;    // response variable\n",
      "    int<lower=1> p_alpha_df;             // prior degrees of freedom for alpha\n",
      "    real p_alpha_loc;                    // prior location for alpha\n",
      "    real<lower=0> p_alpha_scale;         // prior scale for alpha\n",
      "    int<lower=1> p_beta_df;              // prior degrees of freedom for beta\n",
      "    real p_beta_loc;                     // prior location for beta\n",
      "    real<lower=0> p_beta_scale;          // prior scale for beta\n",
      "}\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "parameters {\n",
      "    real alpha;      // intercept\n",
      "    vector[d] beta;  // explanatory variable weights\n",
      "}\n",
      "transformed parameters {\n",
      "    // linear predictor\n",
      "    vector[n] eta;\n",
      "    eta = alpha + X * beta;\n",
      "}\n",
      "model {\n",
      "    alpha ~ student_t(p_alpha_df, p_alpha_loc, p_alpha_scale);\n",
      "    beta ~ student_t(p_beta_df, p_beta_loc, p_beta_scale);\n",
      "    y ~ bernoulli_logit(eta);\n",
      "}\n",
      "generated quantities {\n",
      "    vector[n] log_lik;\n",
      "    for (i in 1:n)\n",
      "        log_lik[i] = bernoulli_logit_lpmf(y[i] | eta[i]);\n",
      "}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "with open('logistic_t.stan') as file:\n",
    "    print(file.read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model = CmdStanModel(stan_file='logistic_t.stan')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set priors and sample from the posterior\n",
    "\n",
    "Here we'll use a Student t prior with 7 degrees of freedom and a scale of 2.5, which, as discussed above, is a reasonable default prior when coefficients should be close to zero but have some chance of being large. CmdStanPy returns the posterior distribution for the parameters describing the uncertainty related to unknown parameter values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "11:15:55 - cmdstanpy - INFO - CmdStan start processing\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e2db6b6dea8c44c9ab4e20ad2b799f61",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "chain 1 |          | 00:00 Status"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cf5582f8da604feab9395c611362d989",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "chain 2 |          | 00:00 Status"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b401afd4bb16492ab445d603bc5aaabc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "chain 3 |          | 00:00 Status"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "caeb8e0bb509471faee90ab25331d5ec",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "chain 4 |          | 00:00 Status"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                                                                                                                                                                                                                                                                                                                "
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "11:15:56 - cmdstanpy - INFO - CmdStan done processing.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "data1 = dict(\n",
    "    n=n,\n",
    "    d=p,\n",
    "    X=X,\n",
    "    y=y,\n",
    "    p_alpha_df=7,\n",
    "    p_alpha_loc=0,\n",
    "    p_alpha_scale=2.5,\n",
    "    p_beta_df=7,\n",
    "    p_beta_loc=0,\n",
    "    p_beta_scale=2.5\n",
    ")\n",
    "fit1 = model.sample(data=data1, seed=74749)\n",
    "idata_t = az.from_cmdstanpy(fit1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inspect the resulting posterior\n",
    "\n",
    "Check n_effs and Rhats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>0.002</td>\n",
       "      <td>0.002</td>\n",
       "      <td>5029.0</td>\n",
       "      <td>3276.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[0]</th>\n",
       "      <td>0.003</td>\n",
       "      <td>0.003</td>\n",
       "      <td>3530.0</td>\n",
       "      <td>3136.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[1]</th>\n",
       "      <td>0.003</td>\n",
       "      <td>0.002</td>\n",
       "      <td>3995.0</td>\n",
       "      <td>3275.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[2]</th>\n",
       "      <td>0.002</td>\n",
       "      <td>0.002</td>\n",
       "      <td>4834.0</td>\n",
       "      <td>3061.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[3]</th>\n",
       "      <td>0.003</td>\n",
       "      <td>0.003</td>\n",
       "      <td>3456.0</td>\n",
       "      <td>2771.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[4]</th>\n",
       "      <td>0.003</td>\n",
       "      <td>0.002</td>\n",
       "      <td>3753.0</td>\n",
       "      <td>3300.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[5]</th>\n",
       "      <td>0.003</td>\n",
       "      <td>0.003</td>\n",
       "      <td>3277.0</td>\n",
       "      <td>2917.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[6]</th>\n",
       "      <td>0.002</td>\n",
       "      <td>0.003</td>\n",
       "      <td>5797.0</td>\n",
       "      <td>3041.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[7]</th>\n",
       "      <td>0.003</td>\n",
       "      <td>0.003</td>\n",
       "      <td>3151.0</td>\n",
       "      <td>3009.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat\n",
       "alpha        0.002    0.002    5029.0    3276.0    1.0\n",
       "beta[0]      0.003    0.003    3530.0    3136.0    1.0\n",
       "beta[1]      0.003    0.002    3995.0    3275.0    1.0\n",
       "beta[2]      0.002    0.002    4834.0    3061.0    1.0\n",
       "beta[3]      0.003    0.003    3456.0    2771.0    1.0\n",
       "beta[4]      0.003    0.002    3753.0    3300.0    1.0\n",
       "beta[5]      0.003    0.003    3277.0    2917.0    1.0\n",
       "beta[6]      0.002    0.003    5797.0    3041.0    1.0\n",
       "beta[7]      0.003    0.003    3151.0    3009.0    1.0"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(idata_t, var_names=['alpha', 'beta'], kind='diagnostics')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "n_effs are high and Rhats<1.1, which is good.\n",
    "\n",
    "Next we check divergences, E-BMFI and treedepth exceedences as explained in [Robust Statistical Workflow with PyStan Case Study](http://mc-stan.org/users/documentation/case-studies/pystan_workflow.html) by Michael Betancourt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "check_div(idata_t)\n",
    "check_energy(idata_t)\n",
    "check_treedepth(idata_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We don't get any warning message, so everything is fine based on these diagnostics and we can proceed with our analysis.\n",
    "\n",
    "Visualise the marginal posterior distributions of each parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAJZCAYAAADlKTK2AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAUTBJREFUeJzt3XtclGX+//H3DAMSooaWh9CwFCZKrcBzmushLXMz7bQd1HJrWw23SNd0i6wsdfPQmtZmmllmv9TSDm7tdtKINA9glqdKykN2MkEUBmUO9+8PvkwRWgoD18j9ej4ePLru+5655jN3MH7muj73fTksy7IEAAAAI5ymAwAAALAzkjEAAACDSMYAAAAMIhkDAAAwiGQMAADAIJIxAAAAg0jGAAAADCIZAwAAMIhkDAAAwCCSMQDGfPfdd3rooYfUp08ftWnTRp06ddKf//xnrVq16oT6WbRokdxut9xut+69995Kx7N582b97W9/U9euXdW2bVv16tVLEydO1P79+4/5nOeee06XXHKJ2rRpo759+2rRokXHfOwPP/yg1NRUDR8+vNIxlr3PtWvX/ubjevXqJbfbrWXLlpXbv2zZsmAfZT9l575///4aPXq0Fi9erMLCwmP2vXbt2uBzAVQdyRgAIz799FNdeeWVWrRokQ4fPqyLL75YrVq10tq1a3X77bfr8ccfP65+9uzZo2nTpsnhcFQpnv/+97+67rrr9L///U9nnHGGevfuLafTqRdeeEFXXHGFdu3aVeE5L7zwgiZNmqTi4mL94Q9/kMfj0UMPPaT58+cf9TUeeugh+f1+Pfjgg1WKNRRiYmI0aNAgDRo0SP3791dKSooiIiL05ptv6v7771f37t31/PPPixXzgOrnMh0AAPs5cuSI/va3v+nAgQPq37+/Jk+erOjoaEmlSdptt92mJ554QqmpqbrooouO2U8gENC4ceMkSVdeeaWWL19eqXh++OEHjRs3Tj6fTw899JCuu+46SZLf79e4ceP0+uuva/To0Vq6dGkw6fP7/Zo9e7bi4uL0+uuvq2HDhtq/f7/69++vf//73xoyZIgiIyODr/HOO+/o3Xff1d///ne1aNGiUnGGUlxcnKZMmVJh/48//qh58+bp+eef1yOPPKLvv/9eY8eONRAhYB+MjAGoce+8846+++471a9fXw8++GAwEZOkdu3a6Y477pAkPfHEE7/Zz/PPP68NGzZozJgxio+Pr3Q8zz33nIqLi9W1a9dgIiZJEREReuCBB1SvXj199tlnysrKCh7bu3ev8vPzdckll6hhw4aSpEaNGumSSy7RwYMHlZubG3xsYWGhJk6cqOTkZN18882VjrMmNG7cWP/4xz+UkZEhSXrmmWe0YcMGw1EBtRvJGIAa99lnn0mSzjvvPNWvX7/C8a5du0qScnJytG/fvqP28dVXX+mxxx5Tx44ddcMNN1QpnnfffVeSNGDAgArH6tatq169ekkqTSLLHDhwQJLUoEGDco8/9dRTJUkejye4b8aMGfrpp580ceJEuVwnx4TEjTfeqLZt20qS5s2bZzgaoHYjGQNQ48oSlbLE5dfi4uIkSZZlaevWrRWOl00fOhwOPfLII1WqFyssLAzWg7Vp0+aojynb/8tYykbifjkC9svtJk2aSJI++eQT/b//9/80ZMiQYHJzsrjiiisklRbs+3w+w9EAtRfJGIAaVzatt2fPnqMe/+X+b775psLxZ555Rps2bdJdd92lM888s0qx7N27N9g+44wzjvqYZs2aVYilUaNGuvDCC/XBBx/oP//5jwoLC7VixQp98MEHcrvdio+Pl9frVUZGhpo1a6Y777yzSnGacN5550kqTZ6//fZbw9EAtdfJMV4OoFbp3LmznnrqKW3ZskVbt27VueeeW+74Sy+9FGz/+hYLX3zxhR5//HFdeOGFGjp0aJVjKSoqCrZPOeWUoz4mJibmqLHce++9Gjp0qO6+++7gvtjYWD388MOSpPnz5+uLL77Q008/HexDkg4fPqw6depUaUQvFO/995SNUEql07JVTXwBHB3JGIAa16VLF3Xo0EHr16/XiBEjNGHCBHXo0EEHDhzQiy++qFdffVWRkZHyer3lEhafz6dx48bJ6XRq0qRJcjrNDu63bdtWK1as0PLly/XDDz+oWbNmGjRokJo1a6bdu3frySef1IABA9SjRw9J0sKFCzVv3jx9//33io6OVp8+fXTfffeVS3qOV7du3XT66acf8/j//ve/cnVrlcFtLYCaQTIGwIiZM2cqLS1NOTk5GjFiRLljw4YNU3Z2tjZv3lyurqxsNG3MmDE6++yzQxJH3bp1g+3i4mLVq1evwmPKkprY2NgKx+Lj45WWllZh//3336/o6Gj94x//kKTgrSJ69+6tjIwM5ebmatasWdq1a5eWLFlywonlX/7yF3Xq1OmYx9etW1flZCw/Pz/YPlZ9H4CqIxkDYESjRo304osvavXq1fr444914MABNWrUSL1791bbtm3VrVs3SVJSUlLwOWVXM65cuVKZmZnl+iur/frggw80ZMgQSaUjUb/nl7fE+Pbbb496V/nvvvuuwmN/y/Lly7VmzRpNmjRJjRo1kiQ9/fTTio+P1+OPPy6Xy6U+ffro0KFDmjt3rlavXh18v+Fky5YtkkoT1qrcOgTAbyMZA2CMw+HQRRddVOHGrrt379a+fft06qmnVqgnk6Ts7Oxj9rlv375j3g7jaGJjY5WQkKBdu3Zp8+bNR03GNm/eLOnngvbfkpeXpylTpqhTp0666qqrJEk//fST9u3bp379+pW7tUVqaqrmzp2rbdu2hWUy9sYbb0gqrfGLiIgwHA1Qe3E1JYCw88wzz0iSrrvuOkVFRQX3v/baa/r888+P+lM2VXj11VcH9x2vPn36SJJWrFhR4VhRUZFWrlwpSbrkkkt+t68pU6aouLhYDz30UHBfWd1bcXFxuceWbVd1KafqsGjRouD94G699VbD0QC1G8kYACN27NhR4epEn8+np556SosXL1ZCQoL++te/huz13nnnHV166aUaNmxYhWPDhg3TKaecotWrV2vJkiXB/WXrSB48eLDc1OmxrF69Wq+99ppGjhypli1bBvc3atRITZs21dq1a7V79+5g36+88oqk4xtxqyn79u3T5MmTNXHiREnS7bffrpSUFMNRAbUb05QAjFi8eLEWL16s8847T02aNFFJSYk2bdqkn376SQkJCZo/f36520FU1aFDh/T111+rpKSkwrEmTZpo8uTJGj16tDIyMvTyyy8rPj5en332mfbs2aPTTjtN06dP/80RrMOHD2vChAlKSkrSn//85wrHR44cqfvvv19XXXWVOnXqpJ07d+rLL79USkqKOnfuHLL3ebzy8/OD63oGAgEVFRVp9+7d2rFjhwKBgGJiYjR69GjdeOONNR4bYDckYwCM6NGjh/bu3autW7dq8+bNioqK0llnnaVbbrlFN910U7n1KmvCZZddphYtWmjOnDnasGGDtm7dqsaNG+vGG2/UyJEjddppp/3m85944gl98803eumll8otEF7muuuuU2RkpObPn69Vq1apXr16uu666/T3v//dyDSlx+MJLqweGRmpunXrqlGjRrrsssvUqVMnXX755Ue9ehRA6DksbiQDAABgDDVjAAAABpGMAQAAGEQyBgAAYBDJGAAAgEEkYwAAAAaRjAEAABhEMgYAAGAQN30FYEuWZUme/aUbMY3Ccn1IAPZAMgbAnrwe+aa0kSS5MnKlqLqGAwJgV0xTAgAAGMRySAAAAAYxMgYAAGAQyRgAAIBBFPADsCXLd0SBtx+WJDn73ieHq47hiADYFSNjAOwp4FNgzVwF1syVAj7T0QCwMUbGANiTM1LOi+8MtgHAFK6mBAAAMIhpSgAAAIOYpgRgS5ZlSV5P6UZkDMshATCGkTEA9uT1yDexlXwTW/2clAGAASRjAAAABlHAD8CWmKYEEC5IxgAAAAximhIAAMAgrqYEYEuWr0SBldMlSc6eo+VwRRmOCIBdkYwBsKeAV4HMmZIkZ4+/SSIZA2AGyRgAe3K65OxyW7ANAKZQwA8AAGAQBfwAAAAGkYwBAAAYRDIGwJaskiJ5M5rKm9FUVkmR6XAA2BjJGAAAgEEU8AOwJcuyJM/+0o2YRiyHBMAYkjEAAACDmKYEAAAwiDsdArAly1eiQNaTkiRnt5EshwTAGJIxAPYU8Crw3hRJkrPrbWI5JACmkIwBsCenS47UG4NtADCFAn4AAACDKOAHAAAwiGQMAADAIJIxALZklRTJ+9BZ8j50FsshATCKqlUA9uUtNh0BAFDAD8CerEBAKvimdKNBczmcTBQAMINkDAAAwCC+CgIAABhEzRgAW7L8XgXWPitJcna6RY6ISMMRAbArkjEA9uQvUeCt+yVJzvY3SiRjAAwhGQNgT44IOdoNDrYBwBQK+AEAAAyigB8AAMAgkjEAAACDSMYA2JJVUiTv5HPlnXwuyyEBMIoCfgD25ckzHQEAUMAPwJ6sQEDa90XpxulJLIcEwBiSMQAAAIPC+qvg2rVr5Xa7NWvWrCr1M2vWLLndbq1duzZEkQEAAIQGNWMAbMnye2VtXCxJclx4HcshATCGZAyAPflL5H9tjCTJ1W4QyyEBMIZkDNVuR65Pu/cEdGYLp1q34lcOYcIRIcc5lwbbAGBKjf/LWFJSosWLF2vVqlXasWOH9u/fr3r16ik1NVUjR47Uueee+7t99OrVS5L06quvatq0aXrvvfd08OBBtWrVSrfeeqsGDBhwzOe+8cYbmjdvnr7++mvVr19fl156qcaMGaPo6OiQxljb5OUFTvg5BwoCenBikTZk+4L72qe6NCGjrk5tcGLlig0bhnV5I05CjshouW5cYDoMAKj5ZKygoECTJk1S+/bt1aNHD9WvX1979uzR+++/r8zMTL3wwgtq167d7/ZTUlKim2++WR6PR1dccYWKi4v11ltvafTo0crPz9eQIUMqPGfRokX68MMP1atXL3Xu3FkffvihFi5cqPz8fE2fPj3kMVYXj6fmL4Dt3jP/hJ/jdEqxsQ7NmBqrlJRI5eR49cDEIg0cXKDACeZ269c0POHXr6qYGEeNvyYAwH5q/NYWJSUlys/PV5MmTcrt//LLL3Xttdfqggsu0LPPPiup9GrKoUOHKi0tTaNGjQo+tlevXtq7d686dOig+fPnKyoqSpL0/fff68orr1RRUZHefffd4GvMmjVLs2fPVr169bRkyRKdffbZkqTDhw9r4MCB2r17t1atWhV8/InEaMJ55+839tonasbUWPXrWye4/d//HdHosYUGIzp+WzY1Mh0CAMAGanzuJyoqqkKSI0mJiYnq1KmT1q9fL6/Xe1x9paenBxMxSWratKmGDh2qkpIS/ec//6nw+KFDhwYTMUmKjo7WgAEDFAgEtGXLlmqJ0e5SUsoXRaemUiSN8GCVeOSd3l7e6e1llXhMhwPAxoxUU2/btk3z5s1Tdna2fvrppwqJTX5+vho3bvybfbhcLl144YUV9rdv316StHXr1grHzjvvvAr7mjZtKkk6ePBgyGOsLiam7Dp0qdyyMTk53nIjY9nZlUtiTbxn1HaWdOCbn9sAYEiNJ2M5OTkaNmyYJOmiiy5Sy5YtFRMTI4fDoXfffVfbt29XSUnJ7/YTFxcn51GWL2nUqHRqqbCw4lRYbGxshX0REaVXUQV+UcQUqhiri4lapg9Xxp3wc9LHHNLESUWyrNIRsexsrx6eXKQO7V2aMbXeCfVF/RZCzhWtiNvfCrYBwJQaT8aeeuoplZSUaNGiRcFRrDKffPLJcfeTn5+vQCBQISHbv7+0nupoiVdNx1ibVOZqxn9Nr6ex4wrL1Yh17RypR6fEKi6OqyNhlsMZIUfziqPrAFDTajwZ2717t0499dQKSU5xcfFRpxaPxefzaePGjUpNTS23f8OGDZJUpdtPhCpGu4uLc2runPrcZwwAgN9Q48MT8fHxKigo0Jdffhnc5/f79c9//lN5eSdWl/TYY4+Vmy78/vvv9fzzzysqKkqXX355WMQIqXUrl3r9IYpEDGHF8vsU2PSKAptekeX3/f4TAKCa1Pi/jjfddJOysrJ0ww036LLLLlNUVJTWrVunH374QR07dtS6deuOq5/TTz89eI+xnj17Bu8zduDAAd13331HvRqypmMEEMb8R+R/+Q5Jkiv5UimCLwsAzKjxkbGePXvq8ccfV4sWLfT6669rxYoVOvvss/Xyyy8rPj7+uPuJiorSs88+q44dO+r111/XK6+8oqZNm2r69OlHveGriRgBhDGHU45WF8vR6mLJQQ0jAHNq/KavoVC2HNL7779vOBIAAICq4esgAACAQSRjAAAABpGMAbAlq8Qj7+MXy/v4xSyHBMCok/LyIWrFAFSdJe374uc2ABhyUhbwA0BVWQG/rF0fS5IcCZ3lcEYYjgiAXZGMAQAAGETNGAAAgEEnZc0YAFSV5ffJ+vwdSZLDfYkc3IEfgCFMUwKwJaukSL6JrSRJroxcOaLqGo4IgF3xVRCAPTmccpzZIdgGAFMYGQMAADCIr4MAAAAGkYwBAAAYRDIGwJYsb7F8T/WT76l+srzFpsMBYGMU8AOwJysga++mYBsATCEZA2BPEXUUcdPCYBsATOFqSgAAAIOoGQMAADCIaUoAtmQF/LK+ypIkOc7uJoczwnBEAOyKaUoAtsRySADCBSNjAOzJ4ZSanvdzGwAMYWQMAADAIL4OAgAAGEQyBgAAYBDJGABbsrzF8j0zSL5nBrEcEgCjKOAHYE9WQNbONcE2AJhCMgbAniLqKOK6p4NtADCFqykBAAAMomYMAADAIKYpAdiSFfDL2pMtSXK0SGU5JADGME0JwJZYDglAuGBkDIBNOaSGZ/3cBgBDGBkDAAAwiAJ+AAAAg0jGAAAADCIZA2BLlvewfAtvlG/hjbK8h02HA8DGKOAHYE+WX9YX7wXbAGAKyRgAe4qIUsSgfwXbAGAKV1MCAAAYRM0YAACAQUxTArAlK+CXfthWutEkmeWQABjDNCUAW2I5JADhgpExADblkOo1/bkNAIYwMgYAAGAQBfwAAAAGkYwBAAAYRDIGwJYs72H5XrpVvpduZTkkAEaRjAGwJ8sva8sKWVtWsBwSAKO4mhKAPUVEyTlgUrANAKZwNSUAAIBBTFMCAAAYxDQlAFuyAgEpf2fpRlxLOZx8NwVgBtOUAGyJ5ZAAhAtGxgDYV3R90xEAACNjAAAAJlEkAQAAYBDJGAAAgEEkYwBsyfIdkW/Z3+Rb9jdZviOmwwFgYyRjAOwp4JO1cYmsjUukgM90NABsjKspAdiTM1LOfhnBNgCYwtWUAAAABjFNCQAAYBDTlABsyQoEpMIfSjdim7AcEgBjSMYA2JOvWL6pF0oqXQ5JLIcEwBCSMQD25eQjEIB5FPADAAAYRJEEAACAQSRjAAAABpGMAbAly3dE/jfGyf/GOJZDAmAUyRgAewr4FFi3QIF1C1gOCYBRJ5yMrV27Vm63W7NmzaqOeEJq1qxZcrvdwZ9p06ZVqb9p06aV6+9kOAcAjsEZKWfP0XL2HM1ySACMCrvrunv16iVJev/990PW56BBgxQfH6/U1NQKxwoLCzVr1iy9/fbb2rdvnxo3bqx+/fopLS1NdeuWv+9Q165dVadOHe3du1fLly8PWXwAap7DFaWIXn83HQYAhF8yVh0GDRqkTp06Vdjv8Xh00003adu2berWrZsuv/xybdu2TfPnz9f69eu1aNEi1alTJ/j4rl27qmvXrlq7di3JGIAatyPXp917AjqzhVOtW9ni4xuwBVv/Nc+bN0/btm3TbbfdpjFjxgT3T5s2TXPnztWCBQt0++23G4wQQHWxLEs6fLB0I7q+HA5HyF8jLy8Qkn4OFAT04MQibcj+ubatfapLEzLq6tQGVSv9bdiQ0mHAtColYxs2bNDMmTO1efNmRUREqEuXLhozZowSEhLKPW7//v2aM2eOVq5cqe+++05169ZVx44dNWrUKCUlJUmSvvnmG/Xu3Tv4HLfbHWynpaVp1KhRKikp0eLFi7Vq1Srt2LFD+/fvV7169ZSamqqRI0fq3HPPPe7YLcvS0qVLFRMTo5EjR5Y7NnLkSC1atEhLly4lGQNOAh5PJe5dXeJR5PTSzxnv6FwpKibEUUnde+aHpB+nU4qNdWjG1FilpEQqJ8erByYWaeDgAgWqmO+tX9MwJDFWVkxM6JNg4GRT6WTsk08+0Zw5c9S9e3cNGTJEX375pd555x1t2LBBS5YsUYsWLSRJu3fv1pAhQ/T999+rW7du6tOnj/bv36+3335bWVlZWrBggc4//3zVr19faWlpeu655yRJw4YNC75Wx44dJUkFBQWaNGmS2rdvrx49eqh+/fras2eP3n//fWVmZuqFF15Qu3btjiv+nTt36scff1S3bt0UE1P+QzgmJkYpKSnKysrSd999p2bNmlX2NAGoAR265J3wc6IjPPr48tJ29155Ouw/HOKoQicQkB7IqKt+fUvLJvr1rSPLkkaPLaxy35U5d6G0ZVMjo68PhINKJ2NZWVl68MEH9ac//Sm476WXXtKECRP0yCOP6KmnnpIkjR07Vvv27dO8efPUvXv34GNHjBihq666Svfdd5/eeOMN1a9fX6NGjQrWYo0aNarCazZo0ECrVq1SkyZNyu3/8ssvde211+qxxx7Ts88+e1zx79q1S5LUsmXLox5v2bKlsrKytHPnTpIxoBY67D9F7d/4RJLks8K/YiMlpfwVn6mpXAEK1BaV/gRq2bKlrr322nL7rr32Wj377LNatWqV8vLy9P3332vjxo266qqryiViknTWWWcFH//FF18Epyt/S1RUVIVETJISExPVqVMnZWVlyev1KjLy9z+kDh06JEmKjY096vGy/YWFVf/mCaB6mZ5qO5ZQjjrl5HiDI2OSlJ3tDUm/4XruADupdDKWkpIip7N84afT6VRKSop27typ7du3a+fOnZJKa8aOdk+ur776Kvjf40nGJGnbtm2aN2+esrOz9dNPP8nrLf+BlJ+fr8aNG1fiHQE4WYVr3dGHK+NC0k/6mEOaOKlIllU6Ipad7dXDk4vUob1LM6bWq1Lf4XruADupdDJ22mmnHXV/o0al8/+HDh1SQUGBJGnVqlVatWrVMfsqLi4+rtfMyckJ1pJddNFFatmypWJiYuRwOPTuu+9q+/btKikpOa6+6tUr/QA71shX2f5jjZwBOLlZvhIF3p0sSXL2GS+HKyrkrxGqKxX/Nb2exo4rLFcj1rVzpB6dEqu4OK6GBE52lU7Gfvrpp6Pu379/v6TSZKcskcnIyNBNN91U2ZcKeuqpp1RSUqJFixapffv25Y598sknJ9RX2RWfZaN3v1a2/1g1ZQBOcgGvAh/9W5Lk7DVGUuiTsVCJi3Nq7pz63GcMqKUq/deck5OjQCBQbqoyEAgoJydHDodD55xzTjAZ27hx43EnY06ns8LUY5ndu3fr1FNPrZCIFRcXa+vWrScUf8uWLdW4cWPl5OTI4/GUu6LS4/EoJydHzZs3p3gfqK2ckXJeNCLYPhm0buVS61amowAQapUe3965c6eWLFlSbt+SJUu0c+dO/eEPf1DDhg3Vrl07nX/++frPf/6jN998s0IfgUBA69atK7evQYMGys/P15EjRyo8Pj4+XgUFBfryyy+D+/x+v/75z38qL+/ECmUdDoeuueYaeTwePfnkk+WOPfnkk/J4PBUuUABQezhcUYq4dIIiLp1QLVOUAHC8Kj0y1q1bNz388MP64IMPlJiYqC+//FIrV65UXFyc7r333uDjpk+frmHDhik9PV3PPfeczj33XEVHR+vbb7/VJ598ory8PH322WfBx3fu3FmbN2/Wrbfeqvbt2ysyMlIdOnRQhw4ddNNNNykrK0s33HCDLrvsMkVFRWndunX64Ycf1LFjxwqJ3e+59dZb9d5772nu3Lnatm2bzj33XG3dulVZWVlq27ZtuXudAQAAVIdKj4xdcMEFWrBggQoLC7Vw4UKtW7dOffr00eLFi4M3fJWkFi1aaPny5RoxYoQ8Ho+WLVuml156Sdu3b1f79u01Y8aMcv2OHDlS1157rb7++mvNmTNHM2fO1McffyxJ6tmzpx5//HG1aNFCr7/+ulasWKGzzz5bL7/8suLj40/4PcTExOiFF17QsGHDlJubq2effVZfffWVhg8frgULFig6OrqypwdAmLMsS5bfW/pjVeIO/gAQIg6rFn8KzZo1S7Nnz9bzzz9/1IXCK2vt2rUaOnRocJkmACcfq6RIvomlBViujFw5ouoajgiAXdnimuihQ4fK7XZr2rRpVepn2rRpcrvdGjp0aIgiAwAAdlerr43u2LGj0tLSgtupqalV6q9r166qU+fnO2CXrZkJ4CQUGSPXPz4PtgHAlFo9TQkAABDubDFNCQAAEK5q9TQlAByL5StRIHOmJMl58Z3cawyAMSRjAOwp4FVg5XRJkrPbSIXzckgAajeSMQD25HTJ2fHmYBsATKGAHwAAwCAK+AEAAAwiGQMAADCIZAyALVklRfJOaC7vhOaySopMhwPAxqhaBWBfAZ/pCACAAn4A9mQFAlLhD6UbsU3kcDJRAMAMkjEAAACD+CoIAABgEDVjAGzJ8pUo8PFcSZKz820shwTAGJIxAPYU8Crwv4mS9H934icZA2AGyRgAe3K65Ljw2mAbAEyhgB8AAMAgCvgBAAAMIhkDAAAwiGQMgC1ZJUXyPpIk7yNJLIcEwCiqVgHY1+GDpiMAAAr4AdiTFQhI+TtLN+JashwSAGNIxgAAAAziqyAAAIBB1IwBsCXL71Vgw0JJkrP9EDkiIg1HBMCuSMYA2JO/RIEV/5AkOS+8TiIZA2AIyRgAe3JEyHHegGAbAEyhgB8AAMAgCvgBAAAMIhkDAAAwiGQMgC1ZJR55H71A3kcvkFXiMR0OABujgB+ATVnSoe9/bgOAIRTwA7AlK+CXfthWutEkWQ4nV1QCMINkDAAAwKBaXTO2bNkyud3u4E96enql+8rMzCzX15AhQ0IYKQAAsCtb1Iz17t1bycnJSkxMDO7bvXu3XnvtNW3ZskVbtmzRjz/+qPj4eL3//vtH7SMhIUFpaWmSpNmzZ9dI3ACqj+X3ytr0iiTJcf5VLIcEwBhbJGN9+vTR4MGDy+3bsGGDZs+erYiICLVq1Uo//fTTb/aRkJCgUaNGSSIZA2oFf4n8y++SJLna/JHlkIBf2ZHr0+49AZ3ZwqnWrWyRLhhj27PboUMHLV68WOecc46io6PVtm1b0yEBqEmOCDmSegfbQG2Qlxeoch8HCgJ6cGKRNmT7gvvap7o0IaOuTm1Q+eqmhg1rdWVUldg2GWvRooVatGhhOgwAhjgio+Uassh0GKglPJ7wuBaue8/8KvfhdEqxsQ7NmBqrlJRI5eR49cDEIg0cXKBAFXK99WsaVjm26hAT4zAdgn2TMQAAQqVDlzzTIYRMICA9kFFX/frWkST161tHliWNHltYpX7D9Rxt2dTIdAi1+2pKAABw4lJSytdQpqZSU1mdGBkDYEtWiUe+J0prxlx3vCdHVIzhiHAyC5cpuFCNPuXkeIMjY5KUne2tcp/hco7CEckYAJuypLyvf24DVRAOdUeS9OHKuCr3kT7mkCZOKpJllY6IZWd79fDkInVo79KMqfUq3W+4nKNwRDIGwJ5c0Yq49fVgG6gNQnHF4r+m19PYcYXlasS6do7Uo1NiFRdHdVN1IBkDYEsOZ4QcCR1NhwGEnbg4p+bOqc99xmoQZxcAAFTQupVLrVuZjsIeSMYA2JLl98na9qYkyZHcX44IPg4BmGHbT5+8vDw9+uijwW2fz6f8/HyNGzcuuG/s2LFq2JCrP4BayX9E/sV/kSS5MnIlkjEAhtj208fj8Wj58uW/uS8tLY1kDKitHE45WnYJtgHAFNsmY82bN9fnn39uOgwAhjgiT5Hrz8t//4EAUM1s8XVw/PjxcrvdSk9Pr3QfmZmZcrvdcrvdIYwMAADYXa0eGUtOTlZaWlpwOzExsdJ9JSQklOsrPj6+SrEBAABIksOyLG49DcB2LG+xfE8PkCS5/rJCjshTDEcEwK5q9cgYAByTFZC+3/JzGwAMYWQMgC1ZAb+sr7IkSY6zu8nhjDAcEQC7IhkDAAAwyBZXUwIAAIQrasYA2JLl98nasVKS5Gjdk+WQABjDNCUAW7JKiuSbWLoKsisjV46ouoYjAmBXfBUEYE8Opxzx5wfbAGAKI2MAAAAG8XUQAADAIJIxAAAAg0jGANiS5S2Wb+4f5Zv7R1neYtPhALAxCvgB2JMVkLV7fbANAKaQjAGwp4g6irj+2WAbAEzhakoAAACDqBkDAAAwiGlKALZkBfyydn0sSXIkdJbDGWE4IgB2xTQlAFtiOSQA4YKRMQA25ZBOT/q5DQCGMDIGAABgEAX8AAAABpGMAQAAGEQyBsCWLG+xfAuulW/BtSyHBMAoCvgB2JMVkJWbGWwDgCkkYwDsKaKOIq5+ItgGAFO4mhIAAMAgasYAAAAMYpoSgC1ZAb+sbz+VJDnOaMdySACMYZoSgC2xHBKAcMHIGACbckinNv+5DQCGMDIGAABgEAX8AAAABpGMAQAAGEQyBsCWLO9h+RbdLN+im2V5D5sOB4CN1epkbNmyZXK73cGf9PT0SveVmZlZrq8hQ4aEMFIANc7yy9r+X1nb/ytZftPRALAxW1xN2bt3byUnJysxMVGSZFmWMjMz9f777ysnJ0fffvutfD6fEhIS1L9/f91yyy2qU6f88igJCQlKS0uTJM2ePbvG3wOAEIuIUsTAacE2AJhSq6+mXLZsmcaPH6/Jkydr8ODBwf1HjhxRu3btFBUVpY4dOyopKUklJSXKysrSzp071bZtWy1cuFCnnHLKUft1u93q2LGjFi5cWFNvBUCI7Mj1afeegM5s4VTrVrb4PgogzNnyk8jpdOquu+7SDTfcoAYNGgT3e71ejRo1SitXrtSiRYt06623GowSwLHk5QVO+DkHCgJ6cGKRNmT7gvvap7o0IaOuTm1QuYqNhg1rdaUHgBpiy2QsMjJSI0aMOOr+22+/XStXrtT69etJxoDf4fGYGVjv3jP/hJ/jdEqxsQ7NmBqrlJRI5WQf0YIZnyp9qJR78GxZlSihXb+m4Qk/JxRiYrhJLVCb2DIZ+y0uV+kpiYhgnTrg93Tokmc6hOMWCEgPZNRVv76l9aB9ewbUK2ugJKnzf9brsD/mhPs09f63bGpk5HUBVA+SsV955ZVXJEkXXXSR4UgAhFpKSmS57cApDXXgQK0tmwVwkiAZ+4UPPvhAixcvVqtWrXTNNdeYDgcIe6am6So7IpWT4w2OjDmi6mpl6kaNHltY6ThMvX8AtQvJ2P/59NNPlZ6ernr16mnmzJmKiuJSd+D3mKpd+nBl3Ak/J33MIU2cVCTLklJTI5Wd7dXDk4vUob1LM6bWq1Qc1G4BCAWSMUmfffaZ/vznP8vpdGrevHnB+5EBCE+VuYrxX9Praey4wnIjYV07R+rRKbGKi+OqSADm2D4Z++yzzzR8+HAFAgHNnz9f7dq1Mx0SgGoQF+fU3Dn1f77PWLMStfx0rPSuZF05Q47IaNMhArApWydjZYmY3+/XM888o/PPP990SACqWetWLrVuJVklXvleWFa6c+BUs0EBsDXbJmObN2/W8OHD5fP5NG/ePF144YWmQwJQkyKi5LzsoWAbAEyxZTJ24MABDR8+XAcPHlT37t21evVqrV69utxj6tWrp5tvvtlMgACqnSMiUhFd/2I6DACwZzJWWFiogoICSdKHH36oDz/8sMJj4uPjScYAAEC1s2Uy1rx5c33++eemwwBgkBUISAXflG40aC6HkysqAZhhi0+f8ePHy+12Kz09vdJ9ZGZmyu12y+12hzAyAMb4iuWb0VG+GR0lX7HpaADYWK0eGUtOTlZaWlpwuyr3D0tISCjXV3x8fJViAxAGIk8xHQEAyGFZFguzAQAAGGKLaUoAAIBwRTIGAABgEMkYAFuyfEfke3W0fK+OluU7YjocADZGMgbAngI+WdmLZGUvkgI+09EAsLFafTUlAByTM1LO3uOCbQAwhaspAQAADGKaEgAAwCCmKQHYkmVZkmd/6UZMIzkcDrMBAbAtkjEA9uT1yDeljSTJlZErRdU1HBAAu2KaEgAAwCAK+AEAAAxiZAwAAMAgkjEAAACDKOAHYEuW74gCbz8sSXL2vU8OVx3DEQGwK0bGANhTwKfAmrkKrJnLckgAjGJkDIA9OSPlvPjOYBsATOFqSgAAAIOYpgQAADCIaUoAtmRZluT1lG5ExrAcEgBjGBkDYE9ej3wTW8k3sdXPSRkAGEAyBgAAYBAF/ABsiWlKAOGCZAwAAMAgpikBAAAM4mpKALZk+UoUWDldkuTsOVoOV5ThiADYFckYAHsKeBXInClJcvb4mySSMQBmkIwBsCenS84utwXbAGAKBfwAAAAGUcAPAABgEMkYAACAQSRjAGzJKimSN6OpvBlNZZUUmQ4HgI2RjAEAABhEAT8AW7IsS/LsL92IacRySACMqdUjY8uWLZPb7Q7+pKenV7qvzMzMcn0NGTIkhJECqGkOh0OOuqeV/pCIATDIFjfX6d27t5KTk5WYmBjc98EHH+jVV1/Vtm3b9NNPP8nr9apZs2ZKSUnRbbfdprPOOqtcHwkJCUpLS5MkzZ49u0bjBwAAtZctkrE+ffpo8ODB5fZlZmZq06ZNateunRo3biyXy6WvvvpKr776qt544w09/fTT6tKlS/DxCQkJGjVqlCSSMaA2sHwlCmQ9KUlydhsZkuWQduT6tHtPQGe2cKp1K1t8vAIIAdt+WowdO1YZGRkV9q9Zs0Y333yzpk2bpldeecVAZABqRMCrwHtTJEkFybdKkZX/ODxQENCDE4u0IdsX3Nc+1aUJGXV1aoOqVYM0bFirq0kAyMbJWJ06dY66v0uXLmrQoIF2795dwxEBtZ/HE0bXC/kiFHH+DXr19SOafOkheQNHKt2V0ynFxjo0Y2qsUlIilZPj1QMTizRwcIECgaqFuX5Nw6p1EGIxMdTXAaFm22TsWDZu3KiCggKlpqaaDgWodTp0yTMdwq/cG5JeAgHpgYy66te39Etev751ZFnS6LGFVe473M7Zlk2NTIcA1Dq2T8aysrK0ceNGlZSUaNeuXVq5cqXi4uI0fvx406EBOImkpESW205NjTzGIwGgPNsnYx999JHmz58f3E5ISNCMGTPUpk0bg1EBtVO4TblJoRt5ysnxBkfGJCk72xuSfsPxnAEILdsnY/fcc4/uueceFRUVKTc3V0888YSuv/56TZo0SX/84x9NhwfUKuFUb2SVFMk3pY02DpYK//qZFBlT6b7SxxzSxElFsqzSEbHsbK8enlykDu1dmjG1XpXiDKdzBqB62D4ZK1O3bl21a9dOTzzxhK666irdf//9uuiii9SwId9KgVrLWyyHpLg4hxxRlb9q8V/T62nsuMJyNWJdO0fq0SmxiovjakgAv41k7FdcLpc6deqk7du367PPPlOPHj1MhwSgOrhOkevudcF2VcTFOTV3Tn3uMwagUvi0OIoff/xRkhQZSQEuUFs5nE4p7syQ9tm6lUutW4W0SwA2YNvx888+++yo+z/88EO9++67ql+/vi644IKaDQoAANiObUfGrr76aiUlJSkpKUlNmzZVcXGxPv/8c23YsEGRkZGaNGmSYmIqX9ALILxZfq8Ca5+VJDk73SJHBCPhAMywbTJ29913a+3atVq/fr3y8vLkdDrVrFkzXXfddRo2bJhatWKuAajV/CUKvHW/JMnZ/kaJZAyAIbZNxm6//XbdfvvtpsMAYIojQo52g4NtADDFFjVj48ePl9vtVnp6eqX7yMzMlNvtltvtDmFkAExxREbLdc2Tcl3zpByR0abDAWBjtXpkLDk5WWlpacHtxMTESveVkJBQrq/4+PgqxQYAACBJDsuyLNNBAAAA2JUtpikB4NeskiJ5J58r7+RzZZUUmQ4HgI3V6mlKAPhNntAsEg4AVcE0JQBbsgIBad8XpRunJ5XekR8ADCAZAwAAMIivggAAAAZRMwbAliy/V9bGxZIkx4XXsRwSAGNIxgDYk79E/tfGSJJc7QaxHBIAY0jGANiTI0KOcy4NtgHAFAr4AQAADKKAHwAAwCCSMQAAAINIxgDYklXikXd6e3mnt5dV4jEdDgAbo4AfgE1Z0oFvfm4DgCEU8AOwJSvgl/Xtp5Ikxxnt5HByRSUAM0jGAAAADKJmDAAAwCBqxgDYkuX3ydr8miTJ0WagHBF8HAIwg08fAPbkPyL/y3dIklzJl0okYwAM4dMHgD05nHK0ujjYBgBTKOAHAAAwiK+DAAAABpGMAQAAGEQyBsCWrBKPvI9fLO/jF7McEgCjKOAHYFOWtO+Ln9sAYAgF/ABsyQr4Ze36WJLkSOjMckgAjCEZAwAAMIiaMQAAAIOoGQNgS5bfJ+vzdyRJDvclLIcEwBimKQHYklVSJN/EVpIkV0auHFF1DUcEwK74KgjAnhxOOc7sEGwDgCmMjAEAABjE10EAAACDSMYAAAAMqtXJ2LJly+R2u4M/6enple4rMzOzXF9DhgwJYaQAaprlLZbvqX7yPdVPlrfYdDgAbMwWBfy9e/dWcnKyEhMTj/mYgoICDRgwQD/++KO6deumZ555ptzxhIQEpaWlSZJmz55drfECqAFWQNbeTcE2AJhii2SsT58+Gjx48G8+5qGHHlJhYeExjyckJGjUqFGSSMaAcLEj16fdewI6s4VTrVud4MdZRB1F3LQw2AYAU2yRjP2e//3vf1qxYoXuv/9+PfTQQ6bDAWwjL69yI1IHCgJ6cGKRNmT7gvvap7o0IaOuTm1wvNUXTun03qXNAqlhw0qFAgBVZvtkLC8vTw888IAGDhyoHj16mA4HMMrjqdk73XTvmV+p5zmdUmysQzOmxiolJVI5OV49MLFIAwcXKFDJGcf1a8IjG4uJcZgOAUANs30yNmHCBEVEROjee+/VoUOHTIcDGNWhS57pEI5LICA9kFFX/fqWTi/261tHliWNHnvsUoNfc8qvjqevlSSt29cpbN77lk2NTIcAoIbZOhl77bXX9Pbbb+uJJ55QgwYNSMaAk0hKSmS57dTUyGM88uiiIo7oqS63SZI6/2e9DvtjQhYbAJwI2yZjP/zwgx555BENGDBAffr0MR0OEBZqeqquKqNROTne4MiYJGVne0/o+Zbl1OcF7mA7XKYpAdiPbZOx++67Ty6XS/fee6/pUICwUdP1Sh+ujKvU89LHHNLESUWyrNIRsexsrx6eXKQO7V2aMbXeCfS0UpL0rqjVAmCOLZOx5cuXKzMzUzNnzlRDLqECjGnYsHL3nf7X9HoaO66wXI1Y186RenRKrOLiavW9rAHUQrZMxrZu3SpJuvPOO496PCsrS263W+ecc45ee+21mgwNwHGIi3Nq7pz6VbvPGACECVt+el144YXyeDwV9ns8Hr355ptq2rSpunXrpmbNmhmIDsDxat3KpdatKvdcy1ss//M3SJIihr4oR+QpIYwMAI6fLZOx/v37q3///hX2f/PNN3rzzTfVunVrPfLIIwYiA1BjrICsnWuCbQAwxZbJGAAooo4irns62AYAU0jGANiSI8IlR5srTIcBACRjv9S8eXN9/vnnpsMAAAA2YotrwMePHy+326309PRK95GZmSm32y232x3CyACYYgX8Cuxap8CudbICftPhALCxWj0ylpycrLS0tOB2YmJipftKSEgo11d8fHyVYgNgmO+w/PNKpyldGblSVF3DAQGwK4dlWZbpIACgplklHvme6C1Jct3xnhxRrE0JwAySMQAAAINsUTMGAAAQrkjGAAAADCIZA2BLlvewfAtvlG/hjbK8h02HA8DGavXVlABwTJZf1hfvBdsAYArJGAB7iohSxKB/BdsAYApXUwIAABhEzRgAAIBBTFMCsCUr4Jd+2Fa60SRZDmeE2YAA2BbTlABsySopkm9iK0mlyyE5WA4JgCGMjAGwKYdUr+nPbQAwhJExAAAAgyjgBwAAMIhkDAAAwCCSMQC2ZHkPy/fSrfK9dCvLIQEwimQMgD1ZfllbVsjasoLlkAAYxdWUAOwpIkrOAZOCbQAwhaspAQAADGKaEgAAwCCmKQHYkhUISPk7SzfiWsrh5LspADOYpgRgSyyHBCBcMDIGwL6i65uOAAAYGQMAADCJIgkAAACDSMYAAAAMIhkDYEuW74h8y/4m37K/yfIdMR0OABsjGQNgTwGfrI1LZG1cIgV8pqMBYGNcTQnAnpyRcvbLCLYBwBSupgQAADCIaUoAAACDmKYEYEtWICAV/lC6EduE5ZAAGEMyBsCefMXyTb1QUulySGI5JACGkIwBsC8nH4EAzKvVBfzLli3T+PHjg9v9+/fXY489Vqm+cnNz1b9//+B2fHy83n///SrHCAAA7M0WXwt79+6t5ORkJSYmBvf9OlH7teeff16dOnUKbsfFxSktLU2S9Nxzz1VfsAAAwFZskYz16dNHgwcPPuqxskTt1+Lj48ttN2zYUKNGjZIkLV++PPRBAsBJbkeuT7v3BHRmC6dat7LFPy9ASNj+r+W3EjUAtZflO6LAWxMkSc7LHpTDVcdwRDUvLy8Qkn4OFAT04MQibcj+eSWD9qkuTcioq1MbVP0q1YYNudIVtZvtkzEANhXwKbBugST93534qz8Z83jCq0S3e8/8kPTjdEqxsQ7NmBqrlJRI5eR49cDEIg0cXKBACPK99WsaVr2TahYT4zAdAk5itk/Gtm7dqgMHDsjn86l58+bq0qWL4uLiTIcFoLo5I+XsOTrYrgkduuTVyOvUtEBAeiCjrvr1LU1o+/WtI8uSRo8tDEn/J8N527KpkekQcBKzfTK2cOHCctvR0dG644479Je//MVQRABqgsMVpYhefzcdRq2RklI+oU1NZb1P4HjZNhlr3ry5MjIy1K1bNzVt2lQFBQVas2aNZsyYoenTp+uUU07RkCFDTIcJoBYJt+m2UI445eR4gyNjkpSd7Q1Z3+F23oBQs20y1rFjR3Xs2DG4HR0drSuvvFLnnXeerrrqKs2ePVvXX3+9XC7bniKgVrMsSzp8sHQjur4cjuqv+Qm3uqIPV4amJCN9zCFNnFQkyyodEcvO9urhyUXq0N6lGVPrVbn/cDtvQKiRafxKYmKiUlNTtXr1auXm5srtdpsOCUB18Hrkm1T6923X5ZBCdZXiv6bX09hxheVqxLp2jtSjU2IVF8eVkMDvIRk7irIC/uLiYsORAED4i4tzau6c+txnDKgk/lp+xe/3a/PmzZKkM844w3A0AKpNZIxcD+wpbbNGZUi0buVS61amowBOPrYdPy5LuH7J7/dr2rRp2rVrlzp16qTGjRsbiAxATXA4HHJERJb+1EC9GAAci22/Dl511VVyu91yu91q0qSJCgoKtG7dOu3cuVNNmzbVI488YjpEAABgA7ZNxoYPH65PPvlEq1evVkFBgSIjI3XmmWdqxIgRuuWWW9SgQQPTIQKoRpavRIF3J0uSnH3Gy+GKMhwRALuybTJ2zz33mA4BgEkBrwIf/VuS5Ow1RhLJGAAzbFEzNn78eLndbqWnp1e6j7LbXLjdbu3duzeE0QEwwhkp50Uj5LxoRI0thwQAR1OrR8aSk5OVlpYW3E5MTKx0X3FxceX6qlev6jcyBGCOwxWliEsnmA4DAOSwLMsyHQQAAIBd1eqRMQA4FsuypICvdMPp4vYWAIyxRc0YAFTg9cj3QAv5HmgheT2mowFgYyRjAAAABlEzBsCWLMuSDh8s3YiuzzQlAGNIxgAAAAximhIAAMAgrqYEYEuWr0SBzJmSJOfFd7IcEgBjSMYA2FPAq8DK6ZIkZ7eRYjkkAKaQjAGwJ6dLzo43B9sAYAoF/AAAAAZRwA8AAGAQyRgAAIBBJGMAbMkqKZJ3QnN5JzSXVVJkOhwANkbVKgD7KlsoHAAMooAfgC1ZgYBU+EPpRmwTOZxMFAAwg2QMAADAIL4KAgAAGETNGABbsnwlCnw8V5Lk7HwbyyEBMIZkDIA9BbwK/G+iJP3fnfhJxgCYQTIGwJ6cLjkuvDbYBgBTKOAHAAAwiAJ+AAAAg0jGAAAADCIZA2BLVkmRvI8kyftIEsshATCKqlUA9nX4oOkIAIACfgD2ZAUCUv7O0o24liyHBMAYkjEAAACD+CoIAABgEDVjAGzJ8nsV2LBQkuRsP0SOiEjDEQGwK5IxAPbkL1FgxT8kSc4Lr5NIxgAYQjIGwJ4cEXKcNyDYBgBTKOAHAAAwiAJ+AAAAg2p1MrZs2TK53e7gT3p6eqX7yszMLNfXkCFDQhgpAACwK1vUjPXu3VvJyclKTEyscGz//v2aM2eOVq1ape+++04xMTFq2bKlBg4cqBtuuCH4uISEBKWlpUmSZs+eXWOxA6geVolHvn91lSS57lotR1SM4YgA2JUtkrE+ffpo8ODBFfZv27ZNw4cP18GDB9WjRw/169dPHo9Hubm5WrlyZYVkbNSoUZJIxoDawZIOff9z2+Z25Pq0e09AZ7ZwqnUrW/zTAIQN2/7FFRYWauTIkZKkV155Reecc0654z6fz0RYAGqKK1quke8G27VBXl7ghJ9zoCCgBycWaUP2z5957VNdmpBRV6c2OPFKloYNa3X1C1AtbJuMvfjii/r222/1yCOPVEjEJMnlsu2pAWzB4YyQmrUJWX8ej/nRte4980/4OU6nFBvr0IypsUpJiVROjlcPTCzSwMEFCpx4bqf1axqe+JOqKCbGUeOvCYSSbTOON998Uw6HQ/369dNXX32ljz76SIcPH9bZZ5+t7t27KyoqynSIAE4iHbrkmQ6hUgIB6YGMuurXt44kqV/fOrIsafTYwkr1Z+I8bNnUqMZfEwglWyZjJSUl+uKLL9SwYUMtXLhQs2bNUuAXXwFbtGihJ554Qm6322CUAKqT5ffK2vSKJMlx/lW2Xg4pJaX8e09Nte+5AEywZTJWUFAgv9+vAwcO6Mknn9Tf//53DRw4UD6fTy+99JL+/e9/a8SIEXrrrbdUp04d0+ECqA7+EvmX3yVJcrX5Y5WXQzIxPfdrlR2VysnxBkfGJCk721vpGMLhPAAnG1smY2WjYH6/XzfeeKOGDx8ePHbnnXfq66+/1ltvvaX//ve/GjhwoKkwAVQnR4QcSb2D7aoKh7qlD1fGnfBz0scc0sRJRbKs0hGx7GyvHp5cpA7tXZoxtd4J9xcO5wE42dgyGatX7+cPmF69elU43qtXL7311lvavHkzyRhQSzkio+Uassh0GCFVmSsZ/zW9nsaOKyxXI9a1c6QenRKruDiujARqgi2TsZiYGDVp0kQ//PCD6tevX+F42b4jR47UdGgAUKPi4pyaO6c+9xkDDLLt157OnTtLknbs2FHhWNm++Pj4Go0JAExp3cqlXn+IIhEDDLBtMvanP/1JkjR37lwdPHgwuH/fvn16/vnn5XQ61bdvX1PhAahmVolH3se6yPtYF1klHtPhALAx234FSklJ0S233KJnn31WV1xxhXr27Cmfz6f33ntP+/fv1913362zzjrLdJgAqo0l5X39cxsADLFtMiZJ48aNU1JSkhYtWqTly5fL4XAoOTlZDz74oC655BLT4QGoTq5oRdz6erANAKbYOhmTpMGDBx91EXEAtZvDGSFHQkfTYQCAPWrGxo8fL7fbrfT09Er3kZmZKbfbzV35AQBASNXqkbHk5GSlpaUFtxMTEyvdV0JCQrm+uNISOLlZfp+sbW9KkhzJ/eWIqNUfhwDCmMOyLCpXAdiOVVIk38RWkiRXRq4cUXUNRwTArvgqCMCeHE45WnYJtgHAFEbGAAAADOLrIAAAgEEkYwAAAAaRjAGwJctbLO8TveV9orcsb7HpcADYGAX8AOzJCkjfb/m5DQCGUMAPwJasgF/WV1mSJMfZ3eRwRhiOCIBdkYwBAAAYRM0YAACAQdSMAbAly++TtWOlJMnRuifLIQEwhmlKALbEckgAwgVfBQHYk8MpR/z5wTYAmMLIGAAAgEF8HQQAADCIZAwAAMAgkjEAtmR5i+Wb+0f55v6R5ZAAGEUBPwB7sgKydq8PtgHAFJIxAPYUUUcR1z8bbAOAKVxNCQAAYBA1YwAAAAYxTQnAlqyAX9aujyVJjoTOcjgjDEcEwK6YpgRgSyyHBCBcMDIGwKYc0ulJP7cBwBBGxgAAAAyigB8AAMAgkjEAAACDSMYA2JLlLZZvwbXyLbiW5ZAAGEUBPwB7sgKycjODbQAwhWQMgD1F1FHE1U8E2wBgCldTAgAAGETNGAAAgEFMUwKwJSvgl/Xtp5IkxxntWA4JgDG1emRs2bJlcrvdwZ/09PRK95WZmVmuryFDhoQwUgA1zndY/jmXyT/nMsl32HQ0AGzMFiNjvXv3VnJyshITE4P7evXqpb179/7m8xYtWqT27dtLkhISEpSWliZJmj17dvUFC6CGOKRTm//cBgBDbJGM9enTR4MHDy63b+jQoTp06FCFx+bn52vRokVq0KCB2rZtG9yfkJCgUaNGSSIZA2oDR1SMIkdvMB0GANgjGTuam2+++aj758+fL0m64oorVKcOl7sDqD125Pq0e09AZ7ZwqnUr2378A2GHv8ZfefnllyVJV199teFIANhdXl5obkZ7oCCgBycWaUO2L7ivfapLEzLq6tQGVS8dbtiwVpcfA9WOZOwXcnJylJubqzZt2uicc84xHQ6AamR5D8u/5K/y+y35r3xKckWbDqmC7j3zQ9KP0ynFxjo0Y2qsUlIilZPj1QMTizRwcIECIcj31q9pWPVOQiwmhjpAnDxIxn6hbFTsmmuuMRwJgGpn+WVt/6+ckrr+4Scd9seYjqjaBALSAxl11a9vaelFv751ZFnS6LGFIem/Q5e8kPQTSls2NTIdAnDcSMb+T1FRkd566y2dcsopGjBggOlwAFS3iChFDJymCQ8WyheINB1NtUtJKf8eU1Nr/3sGThYkY//nzTfflMfj0aBBgxQbG2s6HADVzBERKUf7mzR+saXxpoM5hlCOOOXkeIMjY5KUne0NWd/hOE0JnExIxv7PK6+8IonCfcBuwrm26MOVcSHpJ33MIU2cVCTLKh0Ry8726uHJRerQ3qUZU+tVuf9wPofAyYBkTNKOHTu0ceNGnX322cGbvAKo3axAQNr3RenG6UlyOMPvisBQXaX4r+n1NHZcYbkasa6dI/XolFjFxYXf+wbshmRM3M4CsCVfsXyz/yBJcmXkSlF1zcZTjeLinJo7pz73GQPClO3/Gr1er1577TVFRkbqyiuvNB0OgJoUY69ap9atXGrdynQUAH7N9snY+++/r7y8PPXt21eNGnEpNGAXjqi6ihy/1XQYACDbFwswRQkAAEyy/cjY3LlzTYcAAABszBYjY+PHj5fb7VZ6enql+8jMzJTb7Zbb7Q5hZABMsbyH5Vs6Ur6lI2V5D5sOB4CN1eqRseTkZKWlpQW3ExMTK91XQkJCub7i4+OrFBsAwyy/rE+XlbYHTjUbCwBbc1iWZZkOAgBqmuX3KrD2WUmSs9MtckSwPBAAM0jGAAAADLJFzRgAAEC4qtU1YwBwLFYgIBV8U7rRoHlYLocEwB5IxgDYk69YvhkdJdX+5ZAAhDeSMQD2FXmK6QgAgAJ+AAAAkyiSAAAAMIhkDAAAwCCSMQC2ZPmOyPfqaPleHS3Ld8R0OABsjGQMgD0FfLKyF8nKXiQFfKajAWBjXE0JwJ6ckXL2HhdsA4ApXE0JAABgENOUAAAABjFNCcCWLMuSPPtLN2IayeFwmA0IgG2RjAGwJ69HviltJLEcEgCzSMaqgcfj0fbt202HAeA3WN5i+b/zSpIiNn4iB0sjAagm55xzjmJiYo55nAL+apCTk6PU1FTTYQAAgDCQnZ2tlJSUYx4nGasGjIwBAIAyjIwBAACEMW5tAQAAYBDJGAAAgEEkYwAAAAaRjAEAABjEfcZQLbZt26a33npLW7Zs0ZYtW5Sfn6+OHTtq4cKFlerv008/1axZs7Rx40b5fD4lJSXp5ptvVv/+/UMcefgpLCzUrFmz9Pbbb2vfvn1q3Lix+vXrp7S0NNWte/w3KnW73cc8NmjQIE2ZMiUU4RoXit+VkpISPf3003r99df13XffqUGDBurZs6fuuusuNWrUqBqjDy9VPZfLli3T+PHjj3n8+eefV6dOnUIVblh67bXXlJ2drc2bN+uLL76Q1+vV5MmTNXjw4BPqJxAIaNGiRVqyZIl27dqlmJgYde3aVenp6WrRokU1RR8+QnEe165dq6FDhx7zeGX+v4QKyRiqxbvvvqs5c+YoMjJSZ511lvLz8yvd18cff6xbb71VUVFRuvzyy1W3bl29/fbbSk9P1/fff6/hw4eHMPLw4vF4dNNNN2nbtm3q1q2bLr/8cm3btk3z58/X+vXrtWjRItWpU+e4+4uPj9egQYMq7E9OTg5l2MaE4nclEAhoxIgRysrK0gUXXKC+fftq165dWrp0qdasWaMlS5aoYcOGNfBuzArl313v3r2P+jsWHx8fypDD0syZM7V3717FxcWpcePG2rt3b6X6uf/++7V06VIlJiZqyJAh+vHHH/XWW2/po48+0uLFi9WyZcvQBh5mQnUeJaljx47q2LFjhf1GPwctoBp88cUX1ubNm62SkhLrxx9/tJKSkqybbrrphPvxer1Wnz59rDZt2lhbt24N7j948KDVt29f67zzzrO++eabUIYeVmbOnGklJSVZU6dOLbd/6tSpVlJSkvXUU08dd1+V/X9wsgjV78rLL79sJSUlWXfffbcVCASC+1988UUrKSnJysjIqJb4w0mozuUrr7xiJSUlWa+88kp1hhvWPvroo+C5mjNnTqXOx5o1a6ykpCTrxhtvtI4cORLcv2rVKispKckaPnx4SGMOR6E4jx9//LGVlJRkPf7449URYpVQM4ZqkZiYqPPOO0+RkZFV6ufjjz/W7t27NWDAgHLfWurVq6e//vWv8nq9Wr58eVXDDUuWZWnp0qWKiYnRyJEjyx0bOXKkYmJitHTpUkPRhZ9Q/a6UndO777673OLhf/rTn9SiRQu98cYbOnz4cOjfQBix899dqHXt2rXKI4Blv5N33nmnoqKigvt79Oihjh07KisrS99++22VXiPcheI8hjOmKRHW1q1bJ0nq1q1bhWNl+9avX1+jMdWUnTt36scff1S3bt0q3Lk5JiZGKSkpysrK0nfffadmzZodV58HDx7U4sWLlZ+frwYNGiglJeU3a8lOJqH4XTly5Ig2bdqks846q8IHv8PhUNeuXbV48WJt3rxZ7du3D1Hk4SfUf3dbt27VgQMH5PP51Lx5c3Xp0kVxcXGhCdYG1q5dG/yb/7Xu3btr3bp1Wrduna688sqaD+4ktHPnTi1YsEBHjhxRkyZN1KVLFzVp0sRoTCRjCGs7d+6UJCUkJFQ4dvrppysmJka7du2q4ahqRtn7OlYtSMuWLZWVlaWdO3cedzK2fft23X///eX2de/eXf/85z9P+sL0UPyu7N69W4FA4DfPedlr1eZkLNR/d7++cCc6Olp33HGH/vKXv1QpTjvweDzat2+fkpKSFBERUeF42f+j2vo5WB1WrFihFStWBLddLpduuukmjR079qjnuCaQjCGsFRYWSiqdHjma2NhYHTp0qCZDqjFl7ys2Nvaox8v2l52j3zN8+HD17dtXLVu2VGRkpL788ks9+eSTyszM1O23367Fixcb+yAKhVD8roT6nJ+sQvV317x5c2VkZKhbt25q2rSpCgoKtGbNGs2YMUPTp0/XKaecoiFDhoQ09trmeH8na+vnYCg1bNhQo0ePVs+ePRUfH6/i4mJt3LhR06dP14IFC+RwODRu3DgjsZGM4ZimTJmikpKS43780KFDa/0VPZURLufxnnvuKbd94YUXas6cORo2bJjWrVun9957T3379g3568K+fn3VWnR0tK688kqdd955uuqqqzR79mxdf/31crn4pwjVLzExUYmJicHtmJgY9enTR+eff76uuOIKLVy4ULfddpuRWQL+AnBMixcvlsfjOe7H9+vXL+RJxO996yssLFSDBg1C+pqhVtnzWDYqcaxRmLL9x/rGfDycTqeuueYarVu3Tjk5OSd1MhaK35WaOOcng+r+u0tMTFRqaqpWr16t3NzcWlO3WB2O93fyWKOY+H2nn366evfuraVLl2rTpk3q1atXjcdAMoZj2rhxo+kQgsndrl271KZNm3LH9u3bJ4/Ho3bt2hmI7PhV9jyW1YKU1e/8Wtn+qibAZYXUJ5IwhqNQ/K60aNFCTqez2s95uKuJv7uy37vi4uIq9VPbxcTE6PTTT9c333wjv99foZSgrFbsaPV9OH6mfx+5tQXCWocOHSRJWVlZFY6V7St7TG3TsmVLNW7cWDk5ORUSJY/Ho5ycHDVv3vy4i/ePZdOmTZJK63tOZqH4XYmOjla7du309ddfV7ippGVZWr16tWJiYiokKLVNdf/d+f1+bd68WZJ0xhlnVLofu+jYsWPwb/7XPvzwQ0m193OwppR9Dpq6fQbJGMKC1+tVbm6udu/eXW5/ly5d1KJFC61YsULbtm0L7j906JCeeuopRUZG1trLuR0Oh6655hp5PB49+eST5Y49+eST8ng8uvbaa8vtLy4uVm5uboV7Dn3++efyer0VXiMnJ0fz5s1TZGSkLr300tC/iRp0or8rP/74o3JzcytMxZWd0xkzZsiyrOD+l156SXv27NEf//hHRUdHV++bMSxU57Is4folv9+vadOmadeuXerUqZMaN25cbe/jZJOXl6fc3Fzl5eWV21/2Ozlz5sxy9acffPCB1q1bp27dutXqe3CdqGOdx6P9PkrSc889p7Vr16ply5Zq27ZtTYRYgcP65acNECK5ubmaO3euJOnw4cN66623dNppp6l79+7Bx/xyLcRvvvlGvXv3Vnx8vN5///1yfR1rWZa9e/fqnnvuqfXLIV1//fXavn27unXrpnPPPVdbt25VVlaW2rZtqxdeeKFcYlC29tqv1wEdN26cVq1apdTUVDVr1kwul0tffvmlPvroIzkcDt1///26/vrrTbzFkDqR35Vx48Zp+fLlFdajCwQCuu2224LLIXXo0EG7d+/W22+/rfj4eC1dutTWyyGdyLl0u93BnyZNmqigoEDr1q3Tzp071bRpU73wwgu1fl3FpUuXKjs7W5L0xRdfaMuWLUpJSQlOK6ampuqaa66RJM2aNUuzZ89WWlqaRo0aVa6f++67L7gcUo8ePbRv3z69+eabqlu3rl566SWdddZZNfvGalgozmOvXr3kcrnUpk0bNWnSRMXFxdq0aZO2bt2q+vXr65lnnjFW9kLNGKrFTz/9VOEO3b/ed7wLU3fu3FkvvviiHn/8cb355pvBBYvHjBlT6xcKj4mJ0QsvvBBcKHzt2rU6/fTTNXz4cN1xxx3HPULTu3dvHTx4UNu3b9fq1avl9Xp12mmn6fLLL9ewYcPCvu7ueIXid8XpdOrf//63nn76ab322mtasGCBTj31VF199dW66667bJGISaE5l8OHD9cnn3yi1atXq6CgQJGRkTrzzDM1YsQI3XLLLWF/8U0oZGdnV/gszMnJKTflWJZE/JaHHnpISUlJWrJkiZ5//nnFxMTokksuUXp6us4888yQxx1uQnEe//SnPykrK0vr16/XgQMH5HQ6dcYZZ2jYsGEaPny4mjZtWi2xHw9GxgAAAAyiZgwAAMAgkjEAAACDSMYAAAAMIhkDAAAwiGQMAADAIJIxAAAAg0jGAAAADCIZAwAAMIhkDAAAwCCSMQAAAINIxgAAAAwiGQMAADDo/wM3xastRKARrgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x590 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = az.plot_forest(idata_t,\n",
    "               combined=True,\n",
    "               var_names=['alpha', 'beta'],\n",
    "               )\n",
    "ax[0].axvline(0, color=\"C1\", linestyle=\"dotted\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use Pareto smoothed importance sampling leave-one-out cross-validation to estimate the predictive performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Computed from 4000 posterior samples and 392 observations log-likelihood matrix.\n",
       "\n",
       "         Estimate       SE\n",
       "elpd_loo  -182.33    11.99\n",
       "p_loo       10.98        -\n",
       "------\n",
       "\n",
       "Pareto k diagnostic values:\n",
       "                         Count   Pct.\n",
       "(-Inf, 0.70]   (good)      392  100.0%\n",
       "   (0.70, 1]   (bad)         0    0.0%\n",
       "   (1, Inf)   (very bad)    0    0.0%"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loo_t = az.loo(idata_t)\n",
    "loo_t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Alternative horseshoe prior on weights\n",
    "\n",
    "In this example, with $n \\gg p$ the difference is small, and thus we don’t expect much difference with a different prior and horseshoe prior is usually more useful for $n<p$.\n",
    "\n",
    "The global scale parameter for horseshoe prior is chosen as recommended by Juho Piironen and Aki Vehtari (2017). On the Hyperprior Choice for the Global Shrinkage Parameter in the Horseshoe Prior. Journal of Machine Learning Research: Workshop and Conference Proceedings (AISTATS 2017 Proceedings), accepted for publication. [arXiv preprint arXiv:1610.05559](http://arxiv.org/abs/1610.05559)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/**\n",
      " * Logistic regression HS-prior\n",
      " *\n",
      " * Priors:\n",
      " *     weights - hierarchical shrinkage\n",
      " *     intercept - student t\n",
      " */\n",
      "\n",
      "data {\n",
      "    int<lower=0> n;                     // number of data points\n",
      "    int<lower=1> d;                     // explanatory variable dimension\n",
      "    matrix[n, d] X;                     // explanatory variable\n",
      "    array[n] int<lower=0, upper=1> y;   // response variable\n",
      "    int<lower=1> p_alpha_df;            // prior alpha degrees of freedom\n",
      "    real p_alpha_loc;                   // prior alpha location\n",
      "    real<lower=0> p_alpha_scale;        // prior scale alpha\n",
      "    int<lower=1> p_beta_df;             // prior beta degrees of freedom\n",
      "    int<lower=1> p_beta_global_df;      // prior beta global degrees of freedom\n",
      "    real<lower=0> p_beta_global_scale;  // prior beta global scale\n",
      "}\n",
      "\n",
      "parameters {\n",
      "\n",
      "    // intercept\n",
      "    real alpha;\n",
      "\n",
      "    // auxiliary variables for the variance parameters\n",
      "    vector[d] z;\n",
      "    vector<lower=0>[d] lambda_r1;\n",
      "    vector<lower=0>[d] lambda_r2;\n",
      "    real<lower=0> tau_r1;\n",
      "    real<lower=0> tau_r2;\n",
      "}\n",
      "\n",
      "transformed parameters {\n",
      "\n",
      "    vector<lower=0>[d] lambda;  // local variance parameter\n",
      "    real<lower=0> tau;          // global variance parameter\n",
      "    vector[d] beta;             // explanatory variable weights\n",
      "    vector[n] eta;              // linear predictor\n",
      "\n",
      "    lambda = lambda_r1 .* sqrt(lambda_r2);\n",
      "    tau = tau_r1 * sqrt(tau_r2);\n",
      "    beta = z .* (lambda*tau);\n",
      "    eta = alpha + X * beta;\n",
      "}\n",
      "\n",
      "model {\n",
      "\n",
      "    // student t prior for intercept\n",
      "    alpha ~ student_t(p_alpha_df, p_alpha_loc, p_alpha_scale);\n",
      "\n",
      "    z ~ normal(0.0, 1.0);\n",
      "\n",
      "    // half t priors for lambdas\n",
      "    lambda_r1 ~ normal(0.0, 1.0);\n",
      "    lambda_r2 ~ inv_gamma(0.5*p_beta_df, 0.5*p_beta_df);\n",
      "\n",
      "    // half t priors for tau\n",
      "    tau_r1 ~ normal(0.0, p_beta_global_scale);\n",
      "    tau_r2 ~ inv_gamma(0.5*p_beta_global_df, 0.5*p_beta_global_df);\n",
      "\n",
      "    // observation model\n",
      "    y ~ bernoulli_logit(eta);\n",
      "}\n",
      "\n",
      "generated quantities {\n",
      "    vector[n] log_lik;\n",
      "    for (i in 1:n)\n",
      "        log_lik[i] = bernoulli_logit_lpmf(y[i] | eta[i]);\n",
      "}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "with open('logistic_hs.stan') as file:\n",
    "    print(file.read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model = CmdStanModel(stan_file='logistic_hs.stan')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "11:15:57 - cmdstanpy - INFO - CmdStan start processing\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fe46bdc22fdd438d86e2c794071b8710",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "chain 1 |          | 00:00 Status"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b7aa22b1a494444faa8f8167d004fa9a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "chain 2 |          | 00:00 Status"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9bbecd154f1f463e897f2d5aa5ff1b69",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "chain 3 |          | 00:00 Status"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4d697a0ee5e44c218cc5cb5290be8223",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "chain 4 |          | 00:00 Status"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                                                                                                                                                                                                                                                                                                                "
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "11:16:01 - cmdstanpy - INFO - CmdStan done processing.\n",
      "11:16:01 - cmdstanpy - WARNING - Non-fatal error during sampling:\n",
      "Exception: bernoulli_logit_lpmf: Logit transformed probability parameter[1] is -nan, but must be not nan! (in 'logistic_hs.stan', line 64, column 4 to column 29)\n",
      "Exception: bernoulli_logit_lpmf: Logit transformed probability parameter[1] is -nan, but must be not nan! (in 'logistic_hs.stan', line 64, column 4 to column 29)\n",
      "Exception: bernoulli_logit_lpmf: Logit transformed probability parameter[1] is -nan, but must be not nan! (in 'logistic_hs.stan', line 64, column 4 to column 29)\n",
      "Consider re-running with show_console=True if the above output is unclear!\n",
      "11:16:01 - cmdstanpy - WARNING - Some chains may have failed to converge.\n",
      "\tChain 1 had 6 divergent transitions (0.6%)\n",
      "\tChain 2 had 6 divergent transitions (0.6%)\n",
      "\tChain 3 had 5 divergent transitions (0.5%)\n",
      "\tChain 4 had 6 divergent transitions (0.6%)\n",
      "\tUse the \"diagnose()\" method on the CmdStanMCMC object to see further information.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "p0 = 2 # prior guess for the number of relevant variables\n",
    "tau0 = p0 / (p - p0) * 1 / np.sqrt(n)\n",
    "data2 = dict(\n",
    "    n=n,\n",
    "    d=p,\n",
    "    X=X,\n",
    "    y=y,\n",
    "    p_alpha_df=7,\n",
    "    p_alpha_loc=0,\n",
    "    p_alpha_scale=2.5,\n",
    "    p_beta_df=1,\n",
    "    p_beta_global_df=1,\n",
    "    p_beta_global_scale=tau0\n",
    ")\n",
    "fit2 = model.sample(data=data2, seed=74749)\n",
    "idata_hs = az.from_cmdstanpy(fit2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the horseshoe prior has shrunk the posterior distribution of irrelevant features closer to zero, without affecting the posterior distribution of the relevant features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>-0.973</td>\n",
       "      <td>0.141</td>\n",
       "      <td>-1.230</td>\n",
       "      <td>-0.704</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.003</td>\n",
       "      <td>4363.0</td>\n",
       "      <td>2806.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[0]</th>\n",
       "      <td>0.185</td>\n",
       "      <td>0.174</td>\n",
       "      <td>-0.084</td>\n",
       "      <td>0.510</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.002</td>\n",
       "      <td>1592.0</td>\n",
       "      <td>2761.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[1]</th>\n",
       "      <td>1.138</td>\n",
       "      <td>0.163</td>\n",
       "      <td>0.842</td>\n",
       "      <td>1.444</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.002</td>\n",
       "      <td>4221.0</td>\n",
       "      <td>3439.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[2]</th>\n",
       "      <td>0.014</td>\n",
       "      <td>0.095</td>\n",
       "      <td>-0.186</td>\n",
       "      <td>0.205</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.002</td>\n",
       "      <td>4643.0</td>\n",
       "      <td>3782.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[3]</th>\n",
       "      <td>0.103</td>\n",
       "      <td>0.146</td>\n",
       "      <td>-0.114</td>\n",
       "      <td>0.416</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.002</td>\n",
       "      <td>2547.0</td>\n",
       "      <td>3390.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[4]</th>\n",
       "      <td>-0.021</td>\n",
       "      <td>0.095</td>\n",
       "      <td>-0.238</td>\n",
       "      <td>0.149</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.002</td>\n",
       "      <td>4373.0</td>\n",
       "      <td>3546.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[5]</th>\n",
       "      <td>0.394</td>\n",
       "      <td>0.194</td>\n",
       "      <td>-0.015</td>\n",
       "      <td>0.703</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.003</td>\n",
       "      <td>2170.0</td>\n",
       "      <td>1258.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[6]</th>\n",
       "      <td>0.289</td>\n",
       "      <td>0.162</td>\n",
       "      <td>-0.020</td>\n",
       "      <td>0.549</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.003</td>\n",
       "      <td>1733.0</td>\n",
       "      <td>1142.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta[7]</th>\n",
       "      <td>0.318</td>\n",
       "      <td>0.198</td>\n",
       "      <td>-0.038</td>\n",
       "      <td>0.644</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.003</td>\n",
       "      <td>1778.0</td>\n",
       "      <td>1307.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\n",
       "alpha   -0.973  0.141  -1.230   -0.704      0.002    0.003    4363.0   \n",
       "beta[0]  0.185  0.174  -0.084    0.510      0.004    0.002    1592.0   \n",
       "beta[1]  1.138  0.163   0.842    1.444      0.003    0.002    4221.0   \n",
       "beta[2]  0.014  0.095  -0.186    0.205      0.001    0.002    4643.0   \n",
       "beta[3]  0.103  0.146  -0.114    0.416      0.003    0.002    2547.0   \n",
       "beta[4] -0.021  0.095  -0.238    0.149      0.001    0.002    4373.0   \n",
       "beta[5]  0.394  0.194  -0.015    0.703      0.004    0.003    2170.0   \n",
       "beta[6]  0.289  0.162  -0.020    0.549      0.004    0.003    1733.0   \n",
       "beta[7]  0.318  0.198  -0.038    0.644      0.005    0.003    1778.0   \n",
       "\n",
       "         ess_tail  r_hat  \n",
       "alpha      2806.0    1.0  \n",
       "beta[0]    2761.0    1.0  \n",
       "beta[1]    3439.0    1.0  \n",
       "beta[2]    3782.0    1.0  \n",
       "beta[3]    3390.0    1.0  \n",
       "beta[4]    3546.0    1.0  \n",
       "beta[5]    1258.0    1.0  \n",
       "beta[6]    1142.0    1.0  \n",
       "beta[7]    1307.0    1.0  "
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(idata_hs, var_names=['alpha', 'beta'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAJZCAYAAADlKTK2AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAUM5JREFUeJzt3XtclGX+//H3DAMioAaWh9CwDCbyUIHiIc31kJa5mXbaDmq5ta2GW6RrukWWlrqltqa1mWam2S+1tIObu500Is0DmOUpk/KQlZkgCoMyw9y/P/gyRmgpDPeN3K/n48Gj677vmWs+jMPdZ67rc9+XwzAMQwAAALCE0+oAAAAA7IxkDAAAwEIkYwAAABYiGQMAALAQyRgAAICFSMYAAAAsRDIGAABgIZIxAAAAC5GMAQAAWIhkDIAlfvjhB40fP169evVS69at1aFDB/35z3/WqlWrTqufhQsXyu12y+1266GHHqp0PJs3b9bf/vY3de7cWW3atFGPHj00YcIEHTx48KTPefnll3XllVeqdevW6t27txYuXHjSx+7fv1/JyckaOnRopWMs+z3Xrl37m4/r0aOH3G63li5dWm7/0qVLA32U/ZS993379tXIkSO1aNEiFRQUnLTvtWvXBp4LIDhIxgCY7osvvtB1112nhQsX6ujRo7riiivUsmVLrV27Vvfcc4+eeeaZU+pn7969mjJlihwOR5Xi+e9//6ubb75Z//vf/3TuueeqZ8+ecjqdeuWVV3Tttddq9+7dFZ7zyiuvaOLEiSoqKtIf/vAHeTwejR8/XnPnzj3ha4wfP14lJSV67LHHqhRrMERERGjAgAEaMGCA+vbtq6SkJIWEhOjdd9/VI488oq5du2r+/PlitTzAHC6rAwBgL8eOHdPf/vY3HTp0SH379tWkSZMUHh4uqTRJu/vuu/Xss88qOTlZl19++Un78fv9GjNmjCTpuuuu07JlyyoVz/79+zVmzBj5fD6NHz9eN998sySppKREY8aM0dtvv62RI0dqyZIlgaSvpKREM2fOVHR0tN5++23FxMTo4MGD6tu3r/79739r0KBBCg0NDbzG+++/rw8++EB///vf1bx580rFGUzR0dGaPHlyhf0//fST5syZo/nz5+uJJ57Qjz/+qNGjR1sQIWAvjIwBMNX777+vH374QfXr19djjz0WSMQkqW3btrr33nslSc8+++xv9jN//nxt2LBBo0aNUmxsbKXjefnll1VUVKTOnTsHEjFJCgkJ0aOPPqp69erpyy+/VGZmZuDYvn37lJeXpyuvvFIxMTGSpIYNG+rKK6/U4cOHlZOTE3hsQUGBJkyYoMTERN1xxx2VjtMMjRo10j/+8Q+lp6dLkl588UVt2LDB4qiA2o9kDICpvvzyS0lSq1atVL9+/QrHO3fuLEnKzs7WgQMHTtjHN998o6efflopKSm69dZbqxTPBx98IEnq169fhWORkZHq0aOHpNIkssyhQ4ckSQ0aNCj3+LPOOkuS5PF4AvumTZumn3/+WRMmTJDLdWZMRtx2221q06aNJGnOnDkWRwPUfiRjAExVlqiUJS6/Fh0dLUkyDENbt26tcLxs+tDhcOiJJ56oUr1YQUFBoB6sdevWJ3xM2f5fxlI2EvfLEbBfbjdu3FiS9Pnnn+v//b//p0GDBgWSmzPFtddeK6m0YN/n81kcDVC7kYwBMFXZtN7evXtPePyX+7/77rsKx1988UVt2rRJ999/v84777wqxbJv375A+9xzzz3hY5o2bVohloYNG+qyyy7Txx9/rP/85z8qKCjQ8uXL9fHHH8vtdis2NlZer1fp6elq2rSp7rvvvirFaYVWrVpJKk2ev//+e4ujAWq3M2PMHECt0bFjRz3//PPasmWLtm7dqosvvrjc8ddeey3Q/vUtFnbs2KFnnnlGl112mQYPHlzlWAoLCwPtunXrnvAxERERJ4zloYce0uDBg/XAAw8E9kVFRenxxx+XJM2dO1c7duzQCy+8EOhDko4ePao6depUaUQvGL/77ykboZRKp2WrmvgCODmSMQCm6tSpk9q3b6/169dr2LBhGjdunNq3b69Dhw7p1Vdf1ZtvvqnQ0FB5vd5yCYvP59OYMWPkdDo1ceJEOZ3WDuy3adNGy5cv17Jly7R//341bdpUAwYMUNOmTbVnzx4999xz6tevn7p16yZJWrBggebMmaMff/xR4eHh6tWrlx5++OFySc+p6tKli84555yTHv/f//5Xrm6tMritBWAekjEApps+fbpSU1OVnZ2tYcOGlTs2ZMgQZWVlafPmzeXqyspG00aNGqULLrggKHFERkYG2kVFRapXr16Fx5QlNVFRURWOxcbGKjU1tcL+Rx55ROHh4frHP/4hSYFbRfTs2VPp6enKycnRjBkztHv3bi1evPi0E8u//OUv6tChw0mPr1u3rsrJWF5eXqB9svo+AMFBMgbAdA0bNtSrr76q1atX67PPPtOhQ4fUsGFD9ezZU23atFGXLl0kSQkJCYHnlF3NuHLlSmVkZJTrr6z26+OPP9agQYMklY5E/Z5f3hLj+++/P+Fd5X/44YcKj/0ty5Yt05o1azRx4kQ1bNhQkvTCCy8oNjZWzzzzjFwul3r16qUjR45o9uzZWr16deD3rUm2bNkiqTRhrcqtQwD8PpIxAJZwOBy6/PLLK9zYdc+ePTpw4IDOOuusCvVkkpSVlXXSPg8cOHDS22GcSFRUlOLi4rR7925t3rz5hMnY5s2bJR0vaP8tubm5mjx5sjp06KDrr79ekvTzzz/rwIED6tOnT7lbWyQnJ2v27Nnatm1bjUzG3nnnHUmlNX4hISEWRwPUblxNCaBGefHFFyVJN998s8LCwgL733rrLX311Vcn/CmbKrzhhhsC+05Vr169JEnLly+vcKywsFArV66UJF155ZW/29fkyZNVVFSk8ePHB/aV1b0VFRWVe2zZdlWXcqoOCxcuDNwP7q677rI4GqD2IxkDYLqdO3dWuDrR5/Pp+eef16JFixQXF6e//vWvQXu9999/X1dddZWGDBlS4diQIUNUt25drV69WosXLw7sL1tH8vDhw+WmTk9m9erVeuuttzR8+HC1aNEisL9hw4Zq0qSJ1q5dqz179gT6fuONNySd2oibWQ4cOKBJkyZpwoQJkqR77rlHSUlJFkcF1H5MUwIw3aJFi7Ro0SK1atVKjRs3VnFxsTZt2qSff/5ZcXFxmjt3brnbQVTVkSNH9O2336q4uLjCscaNG2vSpEkaOXKk0tPT9frrrys2NlZffvml9u7dq7PPPltTp079zRGso0ePaty4cUpISNCf//znCseHDx+uRx55RNdff706dOigXbt26euvv1ZSUpI6duwYtN/zVOXl5QXW9fT7/SosLNSePXu0c+dO+f1+RUREaOTIkbrttttMjw2wI5IxAKbr1q2b9u3bp61bt2rz5s0KCwvT+eefrzvvvFO33357ufUqzXD11VerefPmmjVrljZs2KCtW7eqUaNGuu222zR8+HCdffbZv/n8Z599Vt99951ee+21cguEl7n55psVGhqquXPnatWqVapXr55uvvlm/f3vf7dkmtLj8QQWVg8NDVVkZKQaNmyoq6++Wh06dNA111xzwqtHAVQPh8HNZAAAACxDzRgAAICFSMYAAAAsRDIGAABgIZIxAAAAC5GMAQAAWIhkDAAAwEIkYwAAABbipq8AbM8wDMlzsHQjomGNXC8SQO1FMgYAXo98k1tLklzpOVJYpMUBAbATpikBAAAsxHJIAAAAFmJkDAAAwEIkYwAAABaigB+A7Rm+Y/K/97gkydn7YTlcdSyOCICdMDIGAH6f/Gtmy79mtuT3WR0NAJthZAwAnKFyXnFfoA0AZuJqSgAAAAsxTQkAAGAhpikB2J5hGJLXU7oRGsFySABMxcgYAHg98k1oKd+ElseTMgAwCckYAACAhSjgB2B7TFMCsBLJGAAAgIWYpgQAALAQV1MCsD3DVyz/yqmSJGf3kXK4wiyOCICdkIwBgN8rf8Z0SZKz298kkYwBMA/JGAA4XXJ2ujvQBgAzUcAPAABgIQr4AQAALEQyBgAAYCGSMQC2ZxQXypveRN70JjKKC60OB4DNkIwBAABYiAJ+ALZnGIbkOVi6EdGQ5ZAAmIpkDAAAwEJMUwIAAFiIuxsCsD3DVyx/5nOSJGeX4SyHBMBUJGMA4PfK/+FkSZKz891iOSQAZiIZAwCnS47k2wJtADATBfwAAAAWooAfAADAQiRjAAAAFiIZA2B7RnGhvOPPl3f8+SyHBMB0VKoCgCR5i6yOAIBNUcAPwPYMv1/K/650o0EzOZxMGgAwD8kYAACAhfj6BwAAYCFqxgDYnlHilX/tS5IkZ4c75QgJtTgiAHZCMgYAJcXyr3hEkuRsd5tEMgbARCRjAOAIkaPtwEAbAMxEAT8AAICFKOAHAACwEMkYAACAhUjGANieUVwo76SL5Z10McshATAdBfwAIEmeXKsjAGBTFPADsD3D75cO7CjdOCeB5ZAAmIpkDAAAwEI1+uvf2rVr5Xa7NWPGjCr1M2PGDLndbq1duzZIkQEAAAQHNWMAbM8o8crYuEiS5LjsZpZDAmAqkjEAKClWyVujJEmutgNYDgmAqUjGUO125vi0Z69f5zV36sKWfORQAzlC5LjoqkAbAMxk+v8Zi4uLtWjRIq1atUo7d+7UwYMHVa9ePSUnJ2v48OG6+OKLf7ePHj16SJLefPNNTZkyRR9++KEOHz6sli1b6q677lK/fv1O+tx33nlHc+bM0bfffqv69evrqquu0qhRoxQeHh7UGGub3Fz/aT/nUL5fj00o1IYsX2Bfu2SXxqVH6qwGp1+uGBNTo0sccQZzhIbLdds8q8MAYFOmJ2P5+fmaOHGi2rVrp27duql+/frau3evPvroI2VkZOiVV15R27Ztf7ef4uJi3XHHHfJ4PLr22mtVVFSkFStWaOTIkcrLy9OgQYMqPGfhwoX65JNP1KNHD3Xs2FGffPKJFixYoLy8PE2dOjXoMVYXj8f8C2C7ds877ec4nVJUlEPTnopSUlKosrO9enRCofoPzJf/9HM7rV8Tc/pPCoKICIclrwsAsAfTb21RXFysvLw8NW7cuNz+r7/+WjfddJMuvfRSvfTSS5JKr6YcPHiwUlNTNWLEiMBje/TooX379ql9+/aaO3euwsLCJEk//vijrrvuOhUWFuqDDz4IvMaMGTM0c+ZM1atXT4sXL9YFF1wgSTp69Kj69++vPXv2aNWqVYHHn06MVmh1yUHLXvt0TXsqSn161wls//d/xzRydIGFEZ2+LZsaWh0CAKAWM33eJywsrEKSI0nx8fHq0KGD1q9fL6/Xe0p9paWlBRIxSWrSpIkGDx6s4uJi/ec//6nw+MGDBwcSMUkKDw9Xv3795Pf7tWXLlmqJ0e6SksoXQicnUxiNmsco9sg7tZ28U9vJKPZYHQ4Am7Gkmnrbtm2aM2eOsrKy9PPPP1dIbPLy8tSoUaPf7MPlcumyyy6rsL9du3aSpK1bt1Y41qpVqwr7mjRpIkk6fPhw0GOsLlZM17XvVLmlYrKzveVGxrKyKp/EWjVNCTswpEPfHW8DgIlMT8ays7M1ZMgQSdLll1+uFi1aKCIiQg6HQx988IG2b9+u4uLi3+0nOjpazhMsWdKwYemUUkFBxamwqKioCvtCQkqvnPL/oogpWDFWFytqmD5ZGX3az0kbdUQTJhbKMEpHxLKyvHp8UqHat3Np2lP1Trs/ardQbVzhCrlnRaANAGYyPRl7/vnnVVxcrIULFwZGscp8/vnnp9xPXl6e/H5/hYTs4MHSeqoTJV5mx1ibVOZKxn9NrafRYwrK1Yh17hiqJydHKTqaKyNRczicIXI0qzjSDgBmMD0Z27Nnj84666wKSU5RUdEJpxZPxufzaePGjUpOTi63f8OGDZJUpdtPBCtGu4uOdmr2rPrcZwwAgN9g+vBEbGys8vPz9fXXXwf2lZSU6J///Kdyc0+vLunpp58uN134448/av78+QoLC9M111xTI2KEdGFLl3r8IYxEDDWWUeKTf9Mb8m96Q0aJ7/efAABBZPr/HW+//XZlZmbq1ltv1dVXX62wsDCtW7dO+/fvV0pKitatW3dK/ZxzzjmBe4x17949cJ+xQ4cO6eGHHz7h1ZBmxwjgDFFyTCWv3ytJciVeJYXwxQGAeUwfGevevbueeeYZNW/eXG+//baWL1+uCy64QK+//rpiY2NPuZ+wsDC99NJLSklJ0dtvv6033nhDTZo00dSpU094w1crYgRwhnA45Wh5hRwtr5Ac1DMCMJfpN30NhrLlkD766COLIwEAAKgavgICAABYiGQMAADAQiRjAGzPKPbI+8wV8j5zBcshATDdGXnJELViAILLkA7sON4GABOdkQX8ABBMhr9Exu7PJEmOuI5yOEMsjgiAnZCMAQAAWIiaMQAAAAudkTVjABBMRolPxlfvS5Ic7ivl4A78AEzENCUA2zOKC+Wb0FKS5ErPkSMs0uKIANgJX/8AwOGU47z2gTYAmImRMQAAAAvxFRAAAMBCJGMAAAAWIhkDYHuGt0i+5/vI93wfGd4iq8MBYDMU8AOA4Zexb1OgDQBmIhkDgJA6Crl9QaANAGbiakoAAAALUTMGAABgIaYpAdie4S+R8U2mJMlxQRc5nCEWRwTATpimBGB7LIcEwEqMjAGAwyk1aXW8DQAmYmQMAADAQnwFBAAAsBDJGAAAgIVIxgDYnuEtku/FAfK9OIDlkACYjgJ+ADD8MnatCbQBwEwkYwAQUkchN78QaAOAmbiaEgAAwELUjAEAAFiIaUoAtmf4S2TszZIkOZonsxwSAFMxTQnA9lgOCYCVGBkDADmkmPOPtwHARIyMAQAAWIgCfgAAAAuRjAEAAFiIZAyA7Rneo/ItuE2+BbfJ8B61OhwANkMBPwAYJTJ2fBhoA4CZSMYAICRMIQP+FWgDgJm4mhIAAMBC1IwBAABYiGlKALZn+Euk/dtKNxonshwSAFMxTQnA9lgOCYCVGBkDADmkek2OtwHARIyMAQAAWIgCfgAAAAuRjAEAAFiIZAyA7Rneo/K9dpd8r93FckgATEcyBgBGiYwty2VsWc5ySABMx9WUABASJme/iYE2AJiJqykBAAAsxDQlAACAhZimBGB7ht8v5e0q3YhuIYeT76kAzMM0JQDbYzkkAFZiZAwAJCm8vtURALApRsYAAAAsRGEEAACAhUjGAAAALEQyBsD2DN8x+Zb+Tb6lf5PhO2Z1OABshmQMAPw+GRsXy9i4WPL7rI4GgM1wNSUAOEPl7JMeaAOAmbiaEgAAwEJMUwIAAFiIaUoAtmf4/VLB/tKNqMYshwTAVCRjAOArku+pyySVLocklkMCYCKSMQCQJCenQwDWoIAfAADAQhRGAAAAWIhkDAAAwEIkYwBsz/AdU8k7Y1TyzhiWQwJgOpIxAPD75F83T/5181gOCYDpTjsZW7t2rdxut2bMmFEd8QTVjBkz5Ha7Az9TpkypUn9Tpkwp19+Z8B4AOAXOUDm7j5Sz+0iWQwJguhp3LXePHj0kSR999FHQ+hwwYIBiY2OVnJxc4VhBQYFmzJih9957TwcOHFCjRo3Up08fpaamKjKy/L2GOnfurDp16mjfvn1atmxZ0OIDYC2HK0whPf5udRgAbKrGJWPVYcCAAerQoUOF/R6PR7fffru2bdumLl266JprrtG2bds0d+5crV+/XgsXLlSdOnUCj+/cubM6d+6stWvXkowBqJKdOT7t2evXec2durClLU7FAE7C1meAOXPmaNu2bbr77rs1atSowP4pU6Zo9uzZmjdvnu655x4LIwRgBsMwpKOHSzfC68vhcJQ7npvrD9prHcr367EJhdqQdbw2rV2yS+PSI3VWg+CX8cbEUBoM1HRVSsY2bNig6dOna/PmzQoJCVGnTp00atQoxcXFlXvcwYMHNWvWLK1cuVI//PCDIiMjlZKSohEjRighIUGS9N1336lnz56B57jd7kA7NTVVI0aMUHFxsRYtWqRVq1Zp586dOnjwoOrVq6fk5GQNHz5cF1988SnHbhiGlixZooiICA0fPrzcseHDh2vhwoVasmQJyRhwBvJ4TvNe1sUehU4tPed4R+ZIYRHlDnftnhes0OR0SlFRDk17KkpJSaHKzvbq0QmF6j8wX/7g5XwB69fEBL/TKoiIcPz+gwCbqXQy9vnnn2vWrFnq2rWrBg0apK+//lrvv/++NmzYoMWLF6t58+aSpD179mjQoEH68ccf1aVLF/Xq1UsHDx7Ue++9p8zMTM2bN0+XXHKJ6tevr9TUVL388suSpCFDhgReKyUlRZKUn5+viRMnql27durWrZvq16+vvXv36qOPPlJGRoZeeeUVtW3b9pTi37Vrl3766Sd16dJFERHlT7wRERFKSkpSZmamfvjhBzVt2rSybxMAC7TvlHtajw8P8eiza0rbXXvk6mjJ0WqIqpTfLz2aHqk+vUtLIPr0riPDkEaOLqiW1zvd96K6bdnU0OoQgBqn0slYZmamHnvsMf3pT38K7Hvttdc0btw4PfHEE3r++eclSaNHj9aBAwc0Z84cde3aNfDYYcOG6frrr9fDDz+sd955R/Xr19eIESMCtVgjRoyo8JoNGjTQqlWr1Lhx43L7v/76a9100016+umn9dJLL51S/Lt375YktWjR4oTHW7RooczMTO3atYtkDKjljpbUVbt3Ppck+Yzqr95ISip/xWZyMldwAnZW6bNOixYtdNNNN5Xbd9NNN+mll17SqlWrlJubqx9//FEbN27U9ddfXy4Rk6Tzzz8/8PgdO3YEpit/S1hYWIVETJLi4+PVoUMHZWZmyuv1KjT0909sR44ckSRFRUWd8HjZ/oKC6vm2CqD6BHtqLtijS9nZ3sDImCRlZXmD2v8v1bRpSgAVVToZS0pKktNZvjDU6XQqKSlJu3bt0vbt27Vr1y5JpTVjJ7on1zfffBP476kkY5K0bds2zZkzR1lZWfr555/l9ZY/ieXl5alRo0aV+I0A1BbBrkv6ZGV00PpKG3VEEyYWyjBKR8Sysrx6fFKh2rdzadpT9YL2OmWo0QJqvkonY2efffYJ9zdsWFoPcOTIEeXn50uSVq1apVWrVp20r6KiolN6zezs7EAt2eWXX64WLVooIiJCDodDH3zwgbZv367i4uJT6qtevdKT3slGvsr2n2zkDEDtYfiK5f9gkiTJ2WusHK6wcseDeUXiv6bW0+gxBeVqxDp3DNWTk6MUHc2Vj4AdVToZ+/nnn0+4/+DBg5JKk52yRCY9PV233357ZV8q4Pnnn1dxcbEWLlyodu3alTv2+eefn1ZfZVd8lo3e/VrZ/pPVlAGoRfxe+T/9tyTJ2WOUpLDffnwVREc7NXtWfe4zBiCg0meA7Oxs+f3+clOVfr9f2dnZcjgcuuiiiwLJ2MaNG085GXM6nRWmHsvs2bNHZ511VoVErKioSFu3bj2t+Fu0aKFGjRopOztbHo+n3BWVHo9H2dnZatasGcX7gB04Q+W8fFigbYYLW7p0YUtTXgpADVfpMfFdu3Zp8eLF5fYtXrxYu3bt0h/+8AfFxMSobdu2uuSSS/Sf//xH7777boU+/H6/1q1bV25fgwYNlJeXp2PHjlV4fGxsrPLz8/X1118H9pWUlOif//yncnNPr8DW4XDoxhtvlMfj0XPPPVfu2HPPPSePx1PhAgUAtZPDFaaQq8Yp5KpxFaYoAaC6VXpkrEuXLnr88cf18ccfKz4+Xl9//bVWrlyp6OhoPfTQQ4HHTZ06VUOGDFFaWppefvllXXzxxQoPD9f333+vzz//XLm5ufryyy8Dj+/YsaM2b96su+66S+3atVNoaKjat2+v9u3b6/bbb1dmZqZuvfVWXX311QoLC9O6deu0f/9+paSkVEjsfs9dd92lDz/8ULNnz9a2bdt08cUXa+vWrcrMzFSbNm3K3esMAACgOlR6ZOzSSy/VvHnzVFBQoAULFmjdunXq1auXFi1aFLjhqyQ1b95cy5Yt07Bhw+TxeLR06VK99tpr2r59u9q1a6dp06aV63f48OG66aab9O2332rWrFmaPn26PvvsM0lS9+7d9cwzz6h58+Z6++23tXz5cl1wwQV6/fXXFRsbe9q/Q0REhF555RUNGTJEOTk5eumll/TNN99o6NChmjdvnsLDwyv79gA4gxiGIaPEW/pjnObd+wGgihxGLT7zzJgxQzNnztT8+fNPuFB4Za1du1aDBw8OLNME4MxmFBfKN6G0gMuVniNHWKTFEQGwE1tcRz148GC53W5NmTKlSv1MmTJFbrdbgwcPDlJkAADA7mr19dQpKSlKTU0NbCcnJ1epv86dO6tOneN3zS5bMxPAGS40Qq5/fBVoA4CZavU0JQAAQE1ni2lKAACAmqpWT1MCwKkwfMXyZ0yXJDmvuI97jQEwFckYAPi98q+cKklydhmu6lwOCQB+jWQMAJwuOVPuCLQBwEwU8AMAAFiIAn4AAAALkYwBAABYiGQMgO0ZxYXyjmsm77hmMooLrQ4HgM1QqQoAkuT3WR0BAJuigB+A7Rl+v1Swv3QjqrEcTiYNAJiHZAwAAMBCfP0DAACwEDVjAGzP8BXL/9lsSZKz490shwTAVCRjAOD3yv+/CZL0f3fiJxkDYB6SMQBwuuS47KZAGwDMRAE/AACAhSjgBwAAsBDJGAAAgIVIxgDYnlFcKO8TCfI+kcBySABMR6UqAEjS0cNWRwDApijgB2B7ht8v5e0q3YhuwXJIAExFMgYAAGAhvv4BAABYiJoxALZnlHjl37BAkuRsN0iOkFCLIwJgJyRjAFBSLP/yf0iSnJfdLJGMATARyRgAOELkaNUv0AYAM1HADwAAYCEK+AEAACxEMgYAAGAhkjEAtmcUe+R98lJ5n7xURrHH6nAA2AwF/AAgQzry4/E2AJiIAn4Atmf4S6T920o3GifK4eSKSgDmIRkDAACwUK2uGVu6dKncbnfgJy0trdJ9ZWRklOtr0KBBQYwUAADYlS1qxnr27KnExETFx8cH9u3Zs0dvvfWWtmzZoi1btuinn35SbGysPvrooxP2ERcXp9TUVEnSzJkzTYkbgDmMEq+MTW9IkhyXXM9ySABMZYtkrFevXho4cGC5fRs2bNDMmTMVEhKili1b6ueff/7NPuLi4jRixAhJJGNArVNSrJJl90uSXK3/yHJIwP/ZmePTnr1+ndfcqQtb2iJlsIRt39n27dtr0aJFuuiiixQeHq42bdpYHRIAqzhC5EjoGWgDZ7rcXH+Vnn8o36/HJhRqQ5YvsK9dskvj0iN1VoOqVTjFxNTqCqlKsW0y1rx5czVv3tzqMADUAI7QcLkGLbQ6DNQiHo+118Z17Z5Xpec7nVJUlEPTnopSUlKosrO9enRCofoPzJe/anme1q+JqVoHQRIR4bA6hADbJmMAAFSX9p1yrQ6hSvx+6dH0SPXpXUeS1Kd3HRmGNHJ0QZX7rinvzZZNDa0OIYCxQgAAUEFSUvnayeRkaimrCyNjAGzPKPbI92xpzZjr3g/lCIuwOCKc6ayeigvG6FN2tjcwMiZJWVneKvcpWf/e1EQkYwAgQ8r99ngbqCKr65E+WRldpeenjTqiCRMLZRilI2JZWV49PqlQ7du5NO2pelXq2+r3piYiGQMAV7hC7no70AbOdFW9YvFfU+tp9JiCcjVinTuG6snJUYqOpsIp2EjGANiewxkiR1yK1WEANUZ0tFOzZ9XnPmMm4Z0FAAAndGFLly5saXUUtR/JGADbM0p8Mra9K0lyJPaVI4RTIwDz2PaMk5ubqyeffDKw7fP5lJeXpzFjxgT2jR49WjExXPUB1Holx1Sy6C+SJFd6jkQyBsBEtj3jeDweLVu27Df3paamkowBduBwytGiU6ANAGaybTLWrFkzffXVV1aHAaAGcITWlevPy37/gQBQDWzxFXDs2LFyu91KS0urdB8ZGRlyu91yu91BjAwAANhdrR4ZS0xMVGpqamA7Pj6+0n3FxcWV6ys2NrZKsQEAAEiSwzAMbjcNwNYMb5F8L/STJLn+slyO0LoWRwTATmr1yBgAnBLDL/245XgbAEzEyBgA2zP8JTK+yZQkOS7oIoczxOKIANgJyRgAAICFbHE1JQAAQE1FzRgA2zNKfDJ2rpQkOS7sznJIAEzFNCUA2zOKC+WbULoasis9R46wSIsjAmAnfP0DAIdTjthLAm0AMBMjYwAAABbiKyAAAICFSMYAAAAsRDIGwPYMb5F8s/8o3+w/yvAWWR0OAJuhgB8ADL+MPesDbQAwE8kYAITUUcgtLwXaAGAmrqYEAACwEDVjAAAAFmKaEoDtGf4SGbs/kyQ54jrK4QyxOCIAdsI0JQDbYzkkAFZiZAwA5JDOSTjeBgATMTIGAABgIQr4AQAALEQyBgAAYCGSMQC2Z3iL5Jt3k3zzbmI5JACmo4AfAAy/jJyMQBsAzEQyBgAhdRRyw7OBNgCYiaspAQAALETNGAAAgIWYpgRge4a/RMb3X0iSHOe2ZTkkAKZimhKA7bEcEgArMTIGAHJIZzU73gYAEzEyBgAAYCEK+AEAACxEMgYAAGAhkjEAtmd4j8q38A75Ft4hw3vU6nAA2EytTsaWLl0qt9sd+ElLS6t0XxkZGeX6GjRoUBAjBWApo0TG9v/K2P5fySixOhoANmOLqyl79uypxMRExcfHS5IMw1BGRoY++ugjZWdn6/vvv5fP51NcXJz69u2rO++8U3XqlF8SJS4uTqmpqZKkmTNnmv47AKhGIWEK6T8l0AYAM9XqqymXLl2qsWPHatKkSRo4cGBg/7Fjx9S2bVuFhYUpJSVFCQkJKi4uVmZmpnbt2qU2bdpowYIFqlu37gn7dbvdSklJ0YIFC8z6VQCYYGeOT3v2+nVec6cubGmL76oAagBbnm2cTqfuv/9+3XrrrWrQoEFgv9fr1YgRI7Ry5UotXLhQd911l4VRAqhOubn+QPtQvl+PTSjUhixfYF+7ZJfGpUfqrAYVqzliYmp1hQcAk9kyGQsNDdWwYcNOuP+ee+7RypUrtX79epIxwGQej3kD9V275wXaIU6/WjX6RrMejZC7S6KyN5bo0QmF6j8wX35/xeeuXxNjSowREdyAFrADWyZjv8XlKn1LQkJYmw4wW/tOuZa8bqjjqOan9JeyJdc1OerTO1KGIY0cXXDCx5sV55ZNDU15HQDWIhn7lTfeeEOSdPnll1scCQAz+evGyPmLgajk5FDrggFgKyRjv/Dxxx9r0aJFatmypW688UarwwFsx6zpP6n86NbRkgitardRfXofv4o6K8t70ueaGSeA2o9k7P988cUXSktLU7169TR9+nSFhXF5O2A2M2ukPlkZHWinjTqiCRMLZRilI2JZWV49PqlQ7du5NO2pepbGCaD2IxmT9OWXX+rPf/6znE6n5syZE7gfGYDa65dXRP5raj2NHlNQrkasc8dQPTk5StHRXDkJoHrZPhn78ssvNXToUPn9fs2dO1dt27a1OiQAJjsrqlj/7jVGBR392njek2reIoL7jAEwja3PNmWJWElJiV588UVdcsklVocEwApGiYwvlipS0h9umypHmK1PjQBMZtszzubNmzV06FD5fD7NmTNHl112mdUhAbBKSJicV48PtAHATLZMxg4dOqShQ4fq8OHD6tq1q1avXq3Vq1eXe0y9evV0xx13WBMgAFM5QkIV0vkvVocBwKZsmYwVFBQoPz9fkvTJJ5/ok08+qfCY2NhYkjEAAFDtbJmMNWvWTF999ZXVYQCoIQy/X8r/rnSjQTM5nFxBCcA8tjjjjB07Vm63W2lpaZXuIyMjQ263W263O4iRAagRfEXyTUuRb1qK5CuyOhoANlOrR8YSExOVmpoa2K7K/cPi4uLK9RUbG1ul2ADUMKF1rY4AgE05DMMwrA4CAADArmwxTQkAAFBTkYwBAABYiGQMgO0ZvmPyvTlSvjdHyvAdszocADZDMgYAfp+MrIUyshZKfp/V0QCwmVp9NSUAnBJnqJw9xwTaAGAmrqYEAACwENOUAAAAFmKaEoDtGYYheQ6WbkQ0lMPhsDYgALZCMgYAXo98k1tLklzpOVJYpMUBAbATpikBAAAsRAE/AACAhRgZAwAAsBDJGAAAgIUo4Adge4bvmPzvPS5JcvZ+WA5XHYsjAmAnjIwBgN8n/5rZ8q+ZzXJIAEzHyBgAOEPlvOK+QBsAzMTVlAAAABZimhIAAMBCTFMCsD3DMCSvp3QjNILlkACYipExAPB65JvQUr4JLY8nZQBgEpIxAAAAC1HAD8D2mKYEYCWSMQAAAAsxTQkAAGAhrqYEYHuGr1j+lVMlSc7uI+VwhVkcEQA7IRkDAL9X/ozpkiRnt79JIhkDYB6SMQBwuuTsdHegDQBmooAfAADAQhTwAwAAWIhkDAAAwEIkYwBszygulDe9ibzpTWQUF1odDgCbIRkDAACwEAX8AGzPMAzJc7B0I6IhyyEBMFWtHhlbunSp3G534CctLa3SfWVkZJTra9CgQUGMFICVHA6HHJFnl/6QiAEwmS1uqNOzZ08lJiYqPj4+sO/jjz/Wm2++qW3btunnn3+W1+tV06ZNlZSUpLvvvlvnn39+uT7i4uKUmpoqSZo5c6ap8QMAgNrLFslYr169NHDgwHL7MjIytGnTJrVt21aNGjWSy+XSN998ozfffFPvvPOOXnjhBXXq1Cnw+Li4OI0YMUISyRhQ2xi+Yvkzn5MkObsM/93lkHbm+LRnr1/nNXfqwpa2OI0CqEa2PYuMHj1a6enpFfavWbNGd9xxh6ZMmaI33njDgsgAmM7vlf/DyZKk/MS7pNATnxoP5fv12IRCbcjyBfa1S3ZpXHqkzmpQ9aqPmJhaXTkC4CRsm4zVqVPnhPs7deqkBg0aaM+ePSZHBOB0eTxBuv7IF6KQS27Vm28f06SrjsjrP3bChzmdUlSUQ9OeilJSUqiys716dEKh+g/Ml99f9TDWr4mpeifVLCKCmjog2GybjJ3Mxo0blZ+fr+TkZKtDAfA72nfKDWJvD/3uI/x+6dH0SPXpXfplrk/vOjIMaeTogqBEENzfp3ps2dTQ6hCAWsf2yVhmZqY2btyo4uJi7d69WytXrlR0dLTGjh1rdWgAaqCkpNBy28nJoSd5JACcGtsnY59++qnmzp0b2I6Li9O0adPUunVrC6MCcCqCPa13KiNT2dnewMiYJGVleYP2+mfCNCWA4LN9Mvbggw/qwQcfVGFhoXJycvTss8/qlltu0cSJE/XHP/7R6vAA/IZg1S8ZxYXyTW6tjQOlgr9+KYVGnPBxaaOOaMLEQhlG6YhYVpZXj08qVPt2Lk17ql6V46AeC7An2ydjZSIjI9W2bVs9++yzuv766/XII4/o8ssvV0wM31QBW/AWySEpOtohR9iJr2r819R6Gj2moFyNWOeOoXpycpSio7kSEkDlkIz9isvlUocOHbR9+3Z9+eWX6tatm9UhAahurrpyPbAu0D6Z6GinZs+qz33GAAQVZ5ET+OmnnyRJoaEU5gJ24HA6pejzTvnxF7Z06cKW1RgQAFux7bj6l19+ecL9n3zyiT744APVr19fl156qblBAQAA27HtyNgNN9yghIQEJSQkqEmTJioqKtJXX32lDRs2KDQ0VBMnTlRExImLeAHULkaJV/61L0mSnB3ulCOEUXEA5rFtMvbAAw9o7dq1Wr9+vXJzc+V0OtW0aVPdfPPNGjJkiFq2ZA4CsI2SYvlXPCJJcra7TSIZA2Ai2yZj99xzj+655x6rwwBQEzhC5Gg7MNAGADPZomZs7NixcrvdSktLq3QfGRkZcrvdcrvdQYwMQE3gCA2X68bn5LrxOTlCw60OB4DN1OqRscTERKWmpga24+PjK91XXFxcub5iY2OrFBsAAIAkOQzDMKwOAgAAwK5sMU0JAL/FKC6Ud9LF8k66WEZxodXhALCZWj1NCQCnzPP7i4QDQHVgmhKA7Rl+v3RgR+nGOQmld+QHAJOQjAEAAFiIr38AAAAWomYMgO0ZJV4ZGxdJkhyX3cxySABMRTIGACXFKnlrlCTJ1XYAyyEBMBXJGAA4QuS46KpAGwDMRAE/AACAhSjgBwAAsBDJGAAAgIVIxgDYnlHskXdqO3mntpNR7LE6HAA2QwE/AMiQDn13vA0AJqKAH4DtGf4SGd9/IUlynNtWDidXVAIwD8kYAACAhagZAwAAsBA1YwBszyjxydj8liTJ0bq/HCGcGgGYhzMOAJQcU8nr90qSXIlXSSRjAEzEGQcAHE45Wl4RaAOAmSjgBwAAsBBfAQEAACxEMgYAAGAhkjEAtmcUe+R95gp5n7mC5ZAAmI4CfgCQIR3YcbwNACaigB+A7Rn+Ehm7P5MkOeI6shwSAFORjAEAAFiImjEAAAALUTMGwPaMEp+Mr96XJDncV7IcEgBTMU0JwPaM4kL5JrSUJLnSc+QIi7Q4IgB2wtc/AHA45TivfaANAGZiZAwAAMBCfAUEAACwEMkYAACAhWp1MrZ06VK53e7AT1paWqX7ysjIKNfXoEGDghgpACsZ3iL5nu8j3/N9ZHiLrA4HgM3YooC/Z8+eSkxMVHx8/Ekfk5+fr379+umnn35Sly5d9OKLL5Y7HhcXp9TUVEnSzJkzqzVeACYz/DL2bQq0AcBMtkjGevXqpYEDB/7mY8aPH6+CgoKTHo+Li9OIESMkkYwBtU3O7hAdueQlNTrbqXND6lgdDgCbsUUy9nv+97//afny5XrkkUc0fvx4q8MBEES5uScf6TqU79djEwq1IcsnKUWS1C65UOPSI3VWA6diYmp1JQeAGsL2yVhubq4effRR9e/fX926dbM6HKDW8nisuYtO1+55Jz3mdEpRUQ5NeypKSUmhys726tEJheo/MF9+v7R+TUy1xBQR4aiWfgGcmWyfjI0bN04hISF66KGHdOTIEavDAWqt9p1yrQ6hAr9fejQ9Ur17uWR884muvEAyHmqvkQ+WFvFXV8xbNjWsln4BnJlsnYy99dZbeu+99/Tss8+qQYMGJGOADSUlhUq+IpW8fLMkKXn4TosjAmA3tk3G9u/fryeeeEL9+vVTr169rA4HqPWqa8rv9/ze6FZ2tle9uzulJq0kSRs3lgSOWRUzAHuxbTL28MMPy+Vy6aGHHrI6FMAWrKqT+mRl9EmPpY06ogkTC2UYkUq++X1lZXn1+KRCtW/n0rSn6lHbBcAUtkzGli1bpoyMDE2fPl0xMXzzBWqz37oi8l9T62n0mAKNHH38tjadO4bqyclRio7mSkoA5rBlMrZ161ZJ0n333XfC45mZmXK73brooov01ltvmRkaABNFRzs1e1Z97czxac9ev85r7tSFLW15WgRgIVuedS677DJ5PJ4K+z0ej9599101adJEXbp0UdOmTS2IDoDZWp7nVYtVt0o5knHeq3KE1rU6JAA2YstkrG/fvurbt2+F/d99953effddXXjhhXriiScsiAyAJQy/jF1rAm0AMJMtkzEAKCekjkJufiHQBgAzkYwBsD1HiEuO1tdaHQYAmyIZ+4VmzZrpq6++sjoMAABgI7a4dnvs2LFyu91KS0urdB8ZGRlyu91yu91BjAxATWD4S+TfvU7+3etk+Et+/wkAEES1emQsMTFRqampge34+PhK9xUXF1eur9jY2CrFBqAG8R1VyZzSaUpXeo4UFmlxQADsxGEYhmF1EABgJaPYI9+zPSVJrns/lCMswuKIANgJyRgAAICFbFEzBgAAUFORjAEAAFiIZAyA7Rneo/ItuE2+BbfJ8B61OhwANlOrr6YEgFNilMjY8WGgDQBmIhkDgJAwhQz4V6ANAGbiakoAAAALUTMGAABgIaYpAdie4S+R9m8r3WicKIczxNqAANgK05QAbM8oLpRvQktJpcshOVgOCYCJGBkDADmkek2OtwHARIyMAQAAWIgCfgAAAAuRjAEAAFiIZAyA7Rneo/K9dpd8r93FckgATEcyBgBGiYwty2VsWc5ySABMx9WUABASJme/iYE2AJiJqykBAAAsxDQlAACAhZimBGB7ht8v5e0q3YhuIYeT76kAzMM0JQDbYzkkAFZiZAwAJCm8vtURALApRsYAAAAsRGEEAACAhUjGAAAALEQyBsD2DN8x+Zb+Tb6lf5PhO2Z1OABshmQMAPw+GRsXy9i4WPL7rI4GgM1wNSUAOEPl7JMeaAOAmbiaEgAAwEJMUwIAAFiIaUoAtmf4/VLB/tKNqMYshwTAVCRjAOArku+pyySVLocklkMCYCKSMQCQJCenQwDWqNUF/EuXLtXYsWMD23379tXTTz9dqb5ycnLUt2/fwHZsbKw++uijKscIAADszRZfBXv27KnExETFx8cH9v06Ufu1+fPnq0OHDoHt6OhopaamSpJefvnl6gsWAADYii2SsV69emngwIEnPFaWqP1abGxsue2YmBiNGDFCkrRs2bLgBwkANdTOHJ/27PXrvOZOXdjSFv/bAExl+7+q30rUANiD4Tsm/4pxkiTn1Y/J4apjcUSnJzfXXy39Hsr367EJhdqQdXxVgnbJLo1Lj9RZDYJzxWlMDFeuArZPxgCgqMCn0HXzJEnHuqZLYWHWBnSaunbPq5Z+nU4pKsqhaU9FKSkpVNnZXj06oVD9B+bLH6T8b/2amOB0dJoiIhyWvC5wIrZPxrZu3apDhw7J5/OpWbNm6tSpk6Kjo60OC4CJOv3hsO5KGC5JmvOHw/IZRy2OqGbw+6VH0yPVp3fpSGGf3nVkGNLI0QVBe432nXKD1tfp2LKpoSWvC5yI7ZOxBQsWlNsODw/Xvffeq7/85S8WRQTAbD4jTM9/da/VYdRISUnl1+pMTmbtTiDYbJuMNWvWTOnp6erSpYuaNGmi/Px8rVmzRtOmTdPUqVNVt25dDRo0yOowAZjAqqmyYKnO0aXsbG9gZEySsrK8Qe3/TH/vgWCwbTKWkpKilJSUwHZ4eLiuu+46tWrVStdff71mzpypW265RS6Xbd8iwDbq1pV09HDpRnh9ORxnVj3RJyurp7QibdQRTZhYKMMoHRHLyvLq8UmFat/OpWlP1QvKa1C7Bdg4GTuZ+Ph4JScna/Xq1crJyZHb7bY6JADVzeuRb2Lp3/qZuBxSdV2R+K+p9TR6TEG5GrHOHUP15OQoRUdzFSQQLCRjJ1BWwF9UVGRxJABgnehop2bPqs99xoBqxl/Vr5SUlGjz5s2SpHPPPdfiaACYIjRCrkf3lrZZo7KCC1u6dGFLq6MAai/bjjOXJVy/VFJSoilTpmj37t3q0KGDGjVqZEFkAMzmcDjkCAkt/TnD6sUAnPls+xXw+uuvl9vtltvtVuPGjZWfn69169Zp165datKkiZ544gmrQwQAADZg22Rs6NCh+vzzz7V69Wrl5+crNDRU5513noYNG6Y777xTDRo0sDpEACYxfMXyfzBJkuTsNVYO15l1B34AZzbbJmMPPvig1SEAqCn8Xvk//bckydljlCSSMQDmsUXN2NixY+V2u5WWllbpPspuc+F2u7Vv374gRgfAcs5QOS8fJuflwyQnd5gHYK5aPTKWmJio1NTUwHZ8fHyl+4qOji7XV716wbnhIQDrOVxhCrlqnNVhALAph2EYhtVBAAAA2FWtHhkDgFNhGIbk95VuOF3c3gKAqWxRMwYAv8nrke/R5vI92lzyeqyOBoDNkIwBAABYiJoxALZnGIZ09HDpRnh9pikBmIpkDAAAwEJMUwIAAFiIqykB2J7hK5Y/Y7okyXnFfSyHBMBUJGMA4PfKv3KqJMnZZbhYDgmAmUjGAMDpkjPljkAbAMxEAT8AAICFKOAHAACwEMkYAACAhUjGANieUVwo77hm8o5rJqO40OpwANgMlaoAIB1fKBwATEYBPwDbM/x+qWB/6UZUYzmcTBoAMA/JGAAAgIX4+gcAAGAhasYA2J7hK5b/s9mSJGfHu1kOCYCpSMYAwO+V/38TJOn/7sRPMgbAPCRjAOB0yXHZTYE2AJiJAn4AAAALUcAPAABgIZIxAAAAC5GMAbA9o7hQ3icS5H0igeWQAJiOSlUAkKSjh62OAIBNUcAPwPYMv1/K21W6Ed2C5ZAAmIpkDAAAwEJ8/QMAALAQNWMAbM8o8cq/YYEkydlukBwhoRZHBMBOSMYAoKRY/uX/kCQ5L7tZIhkDYCKSMQBwhMjRql+gDQBmooAfAADAQhTwAwAAWKhWJ2NLly6V2+0O/KSlpVW6r4yMjHJ9DRo0KIiRAgAAu7JFzVjPnj2VmJio+Pj4CscOHjyoWbNmadWqVfrhhx8UERGhFi1aqH///rr11lsDj4uLi1NqaqokaebMmabFDqD6GcUe+f7VWZLkun+1HGERFkcEwE5skYz16tVLAwcOrLB/27ZtGjp0qA4fPqxu3bqpT58+8ng8ysnJ0cqVKyskYyNGjJBEMgbUPoZ05Mfj7TPczhyf9uz167zmTl3Y0haneeCMZtu/0oKCAg0fPlyS9MYbb+iiiy4qd9zn81kRFgAruMLlGv5BoG2F3Fx/lfs4lO/XYxMKtSHr+PmrXbJL49IjdVaDylelxMTU6ooWwHK2TcZeffVVff/993riiScqJGKS5HLZ9q0BbMXjMSQ5pQatSncclawYHevaPa/KfTidUlSUQ9OeilJSUqiys716dEKh+g/Ml78Kud76NTFVju1UREQ4THkdoKaxbcbx7rvvyuFwqE+fPvrmm2/06aef6ujRo7rgggvUtWtXhYWFWR0iABO075RrdQhB4/dLj6ZHqk/vOpKkPr3ryDCkkaMLqtSvWe/Rlk0NTXkdoKaxZTJWXFysHTt2KCYmRgsWLNCMGTPk/8XXxubNm+vZZ5+V2+22MEoAZnE5vOrbbLkk6d3v+slnnLl34E9KKh97cvKZ+7sAdmHLZCw/P18lJSU6dOiQnnvuOf39739X//795fP59Nprr+nf//63hg0bphUrVqhOnTpWhwugGq1fEyMVexQ69WFJUvort0gWXE0ZrNGn7GxvYGRMkrKyvFXu06xpSsCubJmMlY2ClZSU6LbbbtPQoUMDx+677z59++23WrFihf773/+qf//+VoUJwAQREQ4ZoSEqSegpSaobGSJHqPm1S5+sjK5yH2mjjmjCxEIZRumIWFaWV49PKlT7di5Ne6pepfullguoXrZMxurVO35S6tGjR4XjPXr00IoVK7R582aSMcAGHKHhcg1aaGkMwbhi8V9T62n0mIJyNWKdO4bqyclRio7mikigprJlMhYREaHGjRtr//79ql+/foXjZfuOHTtmdmgAUGnR0U7NnlWf+4wBZxjbflXq2LGjJGnnzp0VjpXti42NNTUmAAiGC1u61OMPYSRiwBnCtsnYn/70J0nS7Nmzdfjw4cD+AwcOaP78+XI6nerdu7dV4QEwkVHskffpTvI+3UlGscfqcADYjG2/NiUlJenOO+/USy+9pGuvvVbdu3eXz+fThx9+qIMHD+qBBx7Q+eefb3WYAExhSLnfHm8DgIlsm4xJ0pgxY5SQkKCFCxdq2bJlcjgcSkxM1GOPPaYrr7zS6vAAmMUVrpC73g60AcBMtk7GJGngwIEnXEQcgH04nCFyxKVYHQYAm7JFzdjYsWPldruVlpZW6T4yMjLkdru5Kz8AAAiqWj0ylpiYqNTU1MB2fHx8pfuKi4sr1xdXWgK1h1Hik7HtXUmSI7GvHCG1+tQIoIZxGIZBtSoAWzOKC+Wb0FKS5ErPkSMs0uKIANgJX/8AwOGUo0WnQBsAzMTIGAAAgIX4CggAAGAhkjEAAAALkYwBsD3DWyTvsz3lfbanDG+R1eEAsBkK+AHA8Es/bjneBgATUcAPwPYMf4mMbzIlSY4LusjhDLE4IgB2QjIGAABgIWrGAAAALETNGADbM0p8MnaulCQ5LuzOckgATMU0JQDbYzkkAFbi6x8AOJxyxF4SaAOAmRgZAwAAsBBfAQEAACxEMgYAAGAhkjEAtmd4i+Sb/Uf5Zv+R5ZAAmI4CfgAw/DL2rA+0AcBMJGMAEFJHIbe8FGgDgJm4mhIAAMBC1IwBAABYiGlKALZn+Etk7P5MkuSI6yiHM8TiiADYCdOUAGyP5ZAAWImRMQCQQzon4XgbAEzEyBgAAICFKOAHAACwEMkYAACAhUjGANie4S2Sb95N8s27ieWQAJiOAn4AMPwycjICbQAwE8kYAITUUcgNzwbaAGAmrqYEAACwEDVjAAAAFmKaEoDtGf4SGd9/IUlynNuW5ZAAmKpWj4wtXbpUbrc78JOWllbpvjIyMsr1NWjQoCBGCsBSvqMqmXW1SmZdLfmOWh0NAJuxxchYz549lZiYqPj4+MC+Hj16aN++fb/5vIULF6pdu3aSpLi4OKWmpkqSZs6cWX3BArCAQzqr2fE2AJjIFslYr169NHDgwHL7Bg8erCNHjlR4bF5enhYuXKgGDRqoTZs2gf1xcXEaMWKEJJIxoLZxhEUodOQGq8MAYFO2SMZO5I477jjh/rlz50qSrr32WtWpwyXuAGqHnTk+7dnr13nNnbqwpW1P/UCNxF/kr7z++uuSpBtuuMHiSADYQW5u9d5k9lC+X49NKNSGLF9gX7tkl8alR+qsBsEtG46JqdVlyEC1IRn7hezsbOXk5Kh169a66KKLrA4HQDXxeH51e0XfUYW8+VdJUsl1z0uucNNi6do9r1r7dzqlqCiHpj0VpaSkUGVne/XohEL1H5gvf5DzwPVrYoLbYSVERFDzhzMPydgvlI2K3XjjjRZHAqA6te+UW247PMSjz675nySp8x9+1tGSCCvCqhZ+v/RoeqT69C4tu+jTu44MQxo5uiDor/Xr99UKWzY1tDoE4LSRjP2fwsJCrVixQnXr1lW/fv2sDgeAiXz+UI3//NFAu7ZJSir/OyUn177fETiTkYz9n3fffVcej0cDBgxQVFSU1eEAqEYnnk67R5I01txQTBlNys72BkbGJCkry1str1MTpimBMxHJ2P954403JFG4D9hBTaor+mRldLX2nzbqiCZMLJRhlI6IZWV59fikQrVv59K0p+oF9bVq0vsKnElIxiTt3LlTGzdu1AUXXBC4ySsA+zD8funAjtKNcxLkcJp3VWB1X4H4r6n1NHpMQbkasc4dQ/Xk5ChFR3P1I1ATkIyJ21kAtucrkm/mHyRJrvQcKSzS2niCKDraqdmz6nOfMaAGs/1fpNfr1VtvvaXQ0FBdd911VocDwCoRtbve6cKWLl3Y0uooAJyI7ZOxjz76SLm5uerdu7caNuSSaMCOHGGRCh271eowANiU7QsGmKIEAABWsv3I2OzZs60OAQAA2JgtRsbGjh0rt9uttLS0SveRkZEht9stt9sdxMgA1ASG96h8S4bLt2S4DO9Rq8MBYDO1emQsMTFRqampge34+PhK9xUXF1eur9jY2CrFBqAGMUpkfLG0tN3/KWtjAWA7DsMwjN9/GADUXkaJV/61L0mSnB3ulCOE5YIAmIdkDAAAwEK2qBkDAACoqWp1zRgAnArD75fyvyvdaNDM1OWQAIBkDAB8RfJNS5FU+5ZDAlDzkYwBgCSF1rU6AgA2RQE/AACAhSiMAAAAsBDJGAAAgIVIxgDYnuE7Jt+bI+V7c6QM3zGrwwFgMyRjAOD3ychaKCNroeT3WR0NAJvhakoAcIbK2XNMoA0AZuJqSgAAAAsxTQkAAGAhpikB2J5hGJLnYOlGREM5HA5rAwJgKyRjAOD1yDe5tSSWQwJgPpKxauDxeLR9+3arwwBwigxvkUp+8EqSQjZ+LgdLIwEIgosuukgRERG/+zgK+KtBdna2kpOTrQ4DAABYKCsrS0lJSb/7OJKxasDIGAAAYGQMAADgDMCtLQAAACxEMgYAAGAhkjEAAAALkYwBAABYiPuModK2bdumFStWaMuWLdqyZYvy8vKUkpKiBQsWVKq/L774QjNmzNDGjRvl8/mUkJCgO+64Q3379g1y5GeGgoICzZgxQ++9954OHDigRo0aqU+fPkpNTVVk5KnflNTtdp/02IABAzR58uRghFsjBeMzVVxcrBdeeEFvv/22fvjhBzVo0EDdu3fX/fffr4YNG1Zj9DVbVd/bpUuXauzYsSc9Pn/+fHXo0CFY4Z4x3nrrLWVlZWnz5s3asWOHvF6vJk2apIEDB55WP36/XwsXLtTixYu1e/duRUREqHPnzkpLS1Pz5s2rKfqaLxjv79q1azV48OCTHq/MvxfJGCrtgw8+0KxZsxQaGqrzzz9feXl5le7rs88+01133aWwsDBdc801ioyM1Hvvvae0tDT9+OOPGjp0aBAjr/k8Ho9uv/12bdu2TV26dNE111yjbdu2ae7cuVq/fr0WLlyoOnXqnHJ/sbGxGjBgQIX9iYmJwQy7RgnGZ8rv92vYsGHKzMzUpZdeqt69e2v37t1asmSJ1qxZo8WLFysmJsaE36ZmCebfa8+ePU/4OYyNjQ1myGeM6dOna9++fYqOjlajRo20b9++SvXzyCOPaMmSJYqPj9egQYP0008/acWKFfr000+1aNEitWjRIriBnyGC9f5KUkpKilJSUirsr9R51QAqaceOHcbmzZuN4uJi46effjISEhKM22+//bT78Xq9Rq9evYzWrVsbW7duDew/fPiw0bt3b6NVq1bGd999F8zQa7zp06cbCQkJxlNPPVVu/1NPPWUkJCQYzz///Cn3Vdl/lzNZsD5Tr7/+upGQkGA88MADht/vD+x/9dVXjYSEBCM9Pb1a4q/JgvXevvHGG0ZCQoLxxhtvVGe4Z5xPP/008P7NmjWrUu/RmjVrjISEBOO2224zjh07Fti/atUqIyEhwRg6dGhQYz6TBOP9/eyzz4yEhATjmWeeCVpc1Iyh0uLj49WqVSuFhoZWqZ/PPvtMe/bsUb9+/cp9o6hXr57++te/yuv1atmyZVUN94xhGIaWLFmiiIgIDR8+vNyx4cOHKyIiQkuWLLEoujNDsD5TZe/zAw88UG7x8D/96U9q3ry53nnnHR09ejT4v0ANxt9r9ercuXOVRwXLPrf33XefwsLCAvu7deumlJQUZWZm6vvvv6/Sa5ypgvH+VgemKWG5devWSZK6dOlS4VjZvvXr15sak5V27dqln376SV26dKlw5+aIiAglJSUpMzNTP/zwg5o2bXpKfR4+fFiLFi1SXl6eGjRooKSkpN+sJTvTBeMzdezYMW3atEnnn39+hZO3w+FQ586dtWjRIm3evFnt2rULUuQ1X7D/Xrdu3apDhw7J5/OpWbNm6tSpk6Kjo4MTrE2tXbs2cK74ta5du2rdunVat26drrvuOvODq0V27dqlefPm6dixY2rcuLE6deqkxo0bV6ovkjFYbteuXZKkuLi4CsfOOeccRUREaPfu3SZHZZ2y3/VkNR0tWrRQZmamdu3adcrJ2Pbt2/XII4+U29e1a1f985//rJVF6MH4TO3Zs0d+v/83/x3KXstOyViw/15/fcFPeHi47r33Xv3lL3+pUpx25fF4dODAASUkJCgkJKTC8bJ/NzudU6vL8uXLtXz58sC2y+XS7bffrtGjR5/wvf8tJGOwXEFBgaTSaY4TiYqK0pEjR8wMyVJlv2tUVNQJj5ftL3vffs/QoUPVu3dvtWjRQqGhofr666/13HPPKSMjQ/fcc48WLVp02ieOmi4Yn6lg/zvUFsH6e23WrJnS09PVpUsXNWnSRPn5+VqzZo2mTZumqVOnqm7duho0aFBQY7eDU/3c2umcGmwxMTEaOXKkunfvrtjYWBUVFWnjxo2aOnWq5s2bJ4fDoTFjxpxWnyRjNjd58mQVFxef8uMHDx5s26twTldNeW8ffPDBctuXXXaZZs2apSFDhmjdunX68MMP1bt376C/LvBbfn0lWnh4uK677jq1atVK119/vWbOnKlbbrlFLhf/m0LNEh8fr/j4+MB2RESEevXqpUsuuUTXXnutFixYoLvvvvu0Zh34lNvcokWL5PF4Tvnxffr0CXrC8Hvf1AoKCtSgQYOgvqYZKvvelo04nGzEpWz/yb75ngqn06kbb7xR69atU3Z2dq1LxoLxmTLj3+FMVN1/r/Hx8UpOTtbq1auVk5NTq2sbq8Opfm5PNrKJyjvnnHPUs2dPLVmyRJs2bVKPHj1O+bkkYza3ceNGq0MIJHe7d+9W69atyx07cOCAPB6P2rZta0FkVVPZ97aspqOsNufXyvZXNSkuK5I+nYTxTBGMz1Tz5s3ldDqr/d/hTGPG32vZZ7OoqKhK/dhRRESEzjnnHH333XcqKSmpUIJQVit2opo/VF1lP7vc2gKWa9++vSQpMzOzwrGyfWWPsYMWLVqoUaNGys7OrpAoeTweZWdnq1mzZqdcvH8ymzZtklRau1PbBOMzFR4errZt2+rbb7+tcGNIwzC0evVqRUREVEhIarvq/nstKSnR5s2bJUnnnntupfuxs5SUlMC54tc++eQTSfY6p5qp7Lx6urfPIBmDabxer3JycrRnz55y+zt16qTmzZtr+fLl2rZtW2D/kSNH9Pzzzys0NNRWl2A7HA7deOON8ng8eu6558ode+655+TxeHTTTTeV219UVKScnJwK9w766quv5PV6K7xGdna25syZo9DQUF111VXB/yUsdrqfqZ9++kk5OTkVpt7K3udp06bJMIzA/tdee0179+7VH//4R4WHh1fvL1PDBOu9LUu4fqmkpERTpkzR7t271aFDBzVq1Kjafo/aIDc3Vzk5OcrNzS23v+xzO3369HJ1qx9//LHWrVunLl261Mh7bdU0J3t/T/TZlaSXX35Za9euVYsWLdSmTZvTei2H8cszDHAacnJyNHv2bEnS0aNHtWLFCp199tnq2rVr4DG/XPfwu+++U8+ePRUbG6uPPvqoXF8nW15l3759evDBB225HNItt9yi7du3q0uXLrr44ou1detWZWZmqk2bNnrllVfKJQFla6X9em3QMWPGaNWqVUpOTlbTpk3lcrn09ddf69NPP5XD4dAjjzyiW265xYpfsdqdzmdqzJgxWrZsWYU15fx+v+6+++7Ackjt27fXnj179N577yk2NlZLlixhOaRKvrdutzvw07hxY+Xn52vdunXatWuXmjRpoldeecWWayguWbJEWVlZkqQdO3Zoy5YtSkpKCkwrJicn68Ybb5QkzZgxQzNnzlRqaqpGjBhRrp+HH344sBxSt27ddODAAb377ruKjIzUa6+9pvPPP9/cX6yGCMb726NHD7lcLrVu3VqNGzdWUVGRNm3apK1bt6p+/fp68cUXT3uqnpoxVNrPP/9c4U7bv953qotQd+zYUa+++qqeeeYZvfvuu4GFh0eNGmXLhcIjIiL0yiuvBBYKX7t2rc455xwNHTpU99577ymPxvTs2VOHDx/W9u3btXr1anm9Xp199tm65pprNGTIkDOyFu9UBeMz5XQ69e9//1svvPCC3nrrLc2bN09nnXWWbrjhBt1///22TMSk4Ly3Q4cO1eeff67Vq1crPz9foaGhOu+88zRs2DDdeeedZ+RFO8GQlZVV4byanZ1dbsqxLFn4LePHj1dCQoIWL16s+fPnKyIiQldeeaXS0tJ03nnnBT3uM0Uw3t8//elPyszM1Pr163Xo0CE5nU6de+65GjJkiIYOHaomTZqcdlyMjAEAAFiImjEAAAALkYwBAABYiGQMAADAQiRjAAAAFiIZAwAAsBDJGAAAgIVIxgAAACxEMgYAAGAhkjEAAAALkYwBAABYiGQMAADAQiRjAAAAFvr/JR+ZOrdWAkgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 600x590 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = az.plot_forest(idata_hs,\n",
    "               combined=True,\n",
    "               var_names=['alpha', 'beta'],\n",
    "               )\n",
    "ax[0].axvline(0, color=\"C1\", linestyle=\"dotted\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute LOO also for the model with Horseshoe prior. Expected log predictive density is higher, but not significantly. This is not surprising as this is a easy data with $n \\gg p$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Computed from 4000 posterior samples and 392 observations log-likelihood matrix.\n",
       "\n",
       "         Estimate       SE\n",
       "elpd_loo  -181.84    11.10\n",
       "p_loo        9.47        -\n",
       "------\n",
       "\n",
       "Pareto k diagnostic values:\n",
       "                         Count   Pct.\n",
       "(-Inf, 0.70]   (good)      392  100.0%\n",
       "   (0.70, 1]   (bad)         0    0.0%\n",
       "   (1, Inf)   (very bad)    0    0.0%"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loo_hs = az.loo(idata_hs)\n",
    "loo_hs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "elpd_diff: 0.4864 (SE 1.53)\n"
     ]
    }
   ],
   "source": [
    "elpd_diff = loo_hs.loo_i - loo_t.loo_i\n",
    "elpd_diff_se = np.sqrt(np.var(elpd_diff, ddof=1)*n)\n",
    "elpd_diff = np.sum(elpd_diff)\n",
    "print(f'elpd_diff: {elpd_diff:.4} (SE {elpd_diff_se:.3})')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pymc",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
